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EXECUTIVE  SUMMARY 


Atmospheric  ducts  occur  when  the  vertical  refractivity  profile  has  a 
negative  gradient.  Such  ducts  occur  frequently  in  many  parts  of  the  world, 
and  depend  on  weather  conditions.  In  an  atmospheric  duct  environment, 
electromagnetic  energy  may  be  propagated  with  little  attenuation  relative  to 
free  space  over  distances  of  hundreds  of  kilometers,  thereby  greatly 
increasing  the  potential  for  interference  to  communications  and  radar 
systems.  Also,  for  a  radar  system,  the  coverage  can  be  altered. 

The  likelihood  of  experiencing  such  a  field  enhancement  is  greater  at 
higher  radio  propagation  frequencies.  Previously  available  computer  models 
for  determining  ducted  fields  were  limited  in  the  frequencies  and  duct  heights 
they  were  capable  of  considering.  A  mathematical  model,  called  DUCT,  has  thus 
been  developed  to  predict  electromagnetic  field  levels  in  a  duct  environment 
for  ducts  at  any  height  in  the  troposphere,  and  for  propagation  frequencies 
through  SHF. 

The  DUCT  model  is  based  on  a  horizontally  homogeneous  waveguide-mode 
formulation,  which  was  developed  utilizing  Fourier  transform  formalism. 
Numerical  difficulties  encountered  by  previous  investigators  have  been 
overcome  by  the  use  of  a  unique  mathematical  formulation  that  (a)  assures 
linear  independence,  even  in  a  numerical  sense,  of  the  homogeneous  form  of  the 
governing  differential  equation;  and  (b)  provides  flexibility  for  judiciously 
choosing  the  particular  solution  to  the  inhomogeneous  form  of  this 
differential  equation.  In  addition,  criteria  have  been  developed  for 
associating  specific  types  of  modal  field  contributions  with  particular 
portions  of  the  eigenvalue  locus  for  the  atmospheric  waveguide,  thereby 
providing  the  potential  for  increased  computational  efficiency. 

Predictions  of  the  DUCT  model  were  compared  with  measurements  at  beyond- 
line-of-sight  distances  in  both  surface  and  elevated  duct  environments  at 
frequencies  between  65  MHz  and  3.3  GHz.  It  appears  to  be  the  first  model 
capable  of  predictions  that  compare  favorably  (within  a  few  dB)  with 
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measurements  performed  i.i  an  elevated  duct  environment  at  frequencies  as  high 
as  2201  MHz. 

3he  effects  on  the  fields  in  a  duct  environment,  as  a  function  of  duct 
height,  duct  size,  source  height,  observer  height  and  propagation  frequency, 
were  ascertained  by  exercising  the  DUCT  model.  Predictions  of  fields  in  the 
presence  of  more  than  one  duct  in  the  atmosphere  are  also  discussed. 

The  DUCT  program  is  deemed  to  be  a  valid  model  for  predicting  and 
studying  beyond-line-of-sight  fields  in  a  tropospheric  duct  environment,  and 
it  provides  prediction  capabilities  for  certain  cases  that  could  not  be 
adequately  treated  by  previously  existing  models. 
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PREFACE 


The  Electromagnetic  Compatibility  Analysis  Center  (ECAC)  is  a  Department 
of  Defense  facility,  established  to  provide  advice  and  assistance  on 
electromagnetic  compatibility  matters  to  the  Secretary  of  Defense,  the  Joint 
Chiefs,  of  Staff,  the  military  departments  and  other  DoD  components,  the 
center,  located  at  North  Severn,  Annapolis,  Maryland  21402,  is  under  the 
policy  control  of  the  Assistant  Secretary  of  Defense  for  Communication, 
Command,  Control,  and  Intelligence  and  the  Chairman,  Joint  Chiefs  of  Staff,  or 
their  designees,  who  jointly  provide  policy  guidance,  assign  projects,  and 
establish  priorities.  ECAC  functions  under  the  executive  direction  of  the 
Secretary  of  the  Air  Force  and  the  management  and  technical  direction  of  the 
Center  are  provided  by  military  and  civil  service  personnel.  The  technical 
support  function  is  provided  through  an  Air  Force-sponsored  contract  with  the 
IIT  Research  Institute  (IITRI). 

To  the  extent  possible,  all  abbreviations  and  symbols  used  in  this  report 
are  taken  from  American  National  Standard  ANSI  (Y10.19  (1969)  "Letter  Symbols 
for  Units  Used  in  Science  and  Technology"  issued  by  the  American  National 
Standards  Institute,  Inc. 

Users  of  this  report  are  invited  to  submit  comments  that  would  be  useful 
in  revising  or  adding  to  this  material  to  the  Director,  ECAC,  North  Severn, 
Annapolis,  Maryland  21402,  Attention:  XM . 
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Section  1 


SECTION  1 
INTRODUCTION 

BACKGROUND 

Duct  Phenomena 

It  is  known  that  the  index  of  refraction,  n,  of  the  troposphere  varies 
with  altitude.  In  a  "standard"  atmosphere,  n  decreases  with  increasing  height 
above  the  ground,  approaching  unity  as  the  altitude  increases.  At  sea  level, 
the  value  of  n  is  generally  in  the  neighborhood  of  1.000301.  Usually  it  is 
more  convenient  to  use  the  quantity  called  "ref ractivity"  than  it  is  to  use 
the  index  of  refraction.  The  ref ractivity ,  N,  is  related  to  n  through  the 
relationship : 

N  =  ( n-1 )  x  106 

Therefore,  the  refractivity  of  the  troposphere  is  approximately  301  at  sea 
level  and  approaches  0  as  the  altitude  increases. 

For  modeling  purposes,  it  is  often  more  convenient  to  consider  'he  earth 
as  flat  and  to  compensate  for  earth  curvature  through  an  "adjustment"  of  the 
refractivity.  This  adjustment  is  accomplished  by  adding  a  term  to  N  which, 
because  of  Snell's  law,  would  cause  a  ray  to  bend  in  such  a  way  that  its 
height  above  the  "flat  earth"  at  each  point  would  be  the  same  as  that  for  a 
ray  in  an  "unadjusted"  refractivity  environment  over  a  curved  earth.  This  new 
refractivity  is  called  the  "modified  refractivity" ,  M,  and  is  given  as: 

z  6 

M  =  N  +  -  x  10 
a 
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where 

a  =  the  radius  of  the  earth. 

z  =  the  height  above  the  ground . 

For  a  standard  atmosphere,  M  increases  with  increasing  height  as  shown  in  the 
left-hand  portion  of  Figure  la. 

Using  a  geometrical  optics  representation  (ray  tracing),  the  right-hand 
portion  of  Figure  la  illustrates  the  manner  in  which  energy  would  be  radiated 
from  a  source  (e.g.,  a  transmitting  antenna)  in  a  refractivity  environment 
characterized  by  the  left-hand  figure.  In  the  illustration,  r  is  the 
horizontal  distance  along  the  ground.  Ihe  rays  shown  are  bent  in  accordance 
with  Snell's  law,  which  states  that  a  ray  traced  from  a  medium  of  lower 
refractivity  to  a  medium  of  higher  refractivity  will  bend  toward  the  normal  to 
the  interface  of  the  two  media.  Similarly,  a  ray  traced  from  a  medium  of 
higher  refractivity  to  a  medium  of  lower  refractivity  will  be  bent  away  from 
the  normal  to  the  interface.  In  this  case,  the  normal  to  the  "interface"  is 
in  the  z-direction,  and  the  rays  in  Figure  la  are  bent  accordingly. 

Notice  that  no  energy  reaches  point  R,  far  from  the  source,  because  R  is 
located  in  the  flat-earth  representation  at  a  point  that  would  be  beyond  the 
horizon  in  the  curved-earth  representation. 

The  situation  depicted  above  assumed  a  "standard"  atmosphere.  Now 
consider  a  case  for  which  part  of  the  M  versus  z  profile  changes  directions, 
such  as  the  situation  shown  in  Figure  1b.  Such  a  situation  can  be  caused  by 
anomalous  weather  conditions.  In  this  case,  a  portion  of  the  rays  emanating 
from  the  source  will,  in  accordance  with  Snell's  law,  be  bent  in  a  manner  that 
will  confine  them  to  remain  within  a  well-defined  "layer"  of  the  atmosphere. 
This  layer  is  referred  to  as  a  "duct”.  Under  certain  conditions,  waves  can 
propagate  within  a  duct  to  great  distances  (i.e.,  to  beyond-line-of-sight 
distances)  with  little  or  no  attenuation  relative  to  free-space  levels. 
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The  upper  boundary  of  the  duct  is  the  height  at  which  the  modified 
refractivity  gradient  (i.e.,  the  ratio  of  the  change  in  M  to  the  change  in  z) 
changes  from  a  negative  value  to  a  positive  value.  This  is  shown  as  point  A 
in  Figure  1.  The  lower  boundary  of  the  duct  is  determined  by  dropping  a 
vertical  line  from  point  A  in  the  figure.  If  the  vertical  line  intersects  the 
profile,  the  value  of  z  at  that  point  of  intersection  is  the  lower  boundary  of 
the  duct  and  the  duct  is  said  to  be  "elevated."  For  the  case  shown  in  Figure 
1b,  the  duct  will  extend  from  the  height  of  point  B  to  the  height  of  point 
A.  If,  on  the  other  hand,  the  modified  refractivity  profile  were  such  that 
the  vertical  line  reached  the  ground  without  intersecting  the  profile  curve, 
the  duct  would  be  called  a  "surface  duct"  or  a  "ground-based  duct." 

Energy  in  a  duct  may  be  propagated  with  little  attenuation  relative  to 
free-space  over  distances  of  hundreds  of  kilometers  and  may  subsequently 
interfere  with  existing  communications  links.  In  order  to  predict  whether 
such  interference  will  take  place,  information  is  required  about  the 
characteristics  of  any  ducts  that  are  present  and  the  effects  of  these  ducts 
on  the  propagated  fields.  Unfortunately,  it  is  not  possible  to  predict  in 
advance  the  occurrence  of  a  duct  at  a  particular  time  and  in  a  given  region. 
However,  statistical  data  is  available  on  the  occurrence  of  elevated  and 
surface  ducts  in  different  months  and  in  different  regions  of  the  world.1 
Reference  1  also  contains  statistical  data  related  to  the  characteristics  of 
these  ducts.  In  general,  the  percent  of  occurrence  of  ducts  can  vary  from  0- 
60%  depending  on  the  region  and  time  of  year.  Therefore,  it  would  be  expected 
that  meaningful  statistical  data  can  be  obtained  for  the  propagated  field 
strengths  observed  in  a  duct  environment.  Ib  accomplish  the  determination  of 
field  strengths,  the  statistical  duct  occurrence  data  would  be  integrated  with 
a  "deterministic"  model  that  calculates  the  field  strength  when  the  duct 
characteristics  are  known.  The  goal  of  this  study  was  to  develop  such  a 
deterministic  model. 


1  Ortenburger,  L.  N. ,  Iawson,  S.  B.  and  Miller,  G.  K. ,  Radiosonde  Data 

Analysis  Summary  Maps  of  Observed  Data,  GTE  Sylvania,  Inc.,  December  1978. 
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Computational  Capabilities 


The  existence  of  ducts  and  their  effects  on  electromagnetic  wave 

propagation  have  been  recognized  for  a  long  time.  Documentation  of  the 

computations  in  this  area  first  appeared  at  the  end  of  World  War  XI.  A  review 

2 

of  the  work  performed  at  that  time  is  available  in  the  text  edited  by  Kerr. 
That  work,  as  in  most  subsequent  studies,  considers  the  propagation 
environment  as  a  waveguide,  the  total  field  strength  at  the  location  of  a 
receiving  antenna  is  then  the  sum  of  the  field  strengths  of  the  modes  of  the 
waveguide.  The  wavequide  is  modeled  as  flat,  with  earth  curvature  compensated 
by  a  modified  refractivity .  justification  for  this  approximation  is  discussed 
by  Pekeris.3,4  The  medium  is  assumed  to  be  laterally  homogeneous. 

Each  waveguide  mode  corresponds  to  an  "eigenvalue"  of  the  system.  In 
many  works,5,6  these  "eigenvalues"  are  associated  with  takeoff  angles  of  a  ray 
relative  to  the  ground  —  waveguide  modes  exist  for  discrete  values  of 


2 

Kerr,  D.  E. ,  Propagation  of  Short  Radio  Waves,  MIT  Radiation  Laboratory 
Series,  Vol.  13,  McGraw-Hill  Book,  New  York,  NY,  1951. 

3  Pekeris,  C.  L. ,  "Wave  Theoretical  Interpretation  of  Propagation  of 
10-Centimeter  and  3-Centimeter  Waves  in  Low-Level  Ocean  Ducts," 

Proc.  of  the  IRE,  May  1947,  pp.  453-462. 

4 

Pekeris,  C.  L. , "Accuracy  of  the  Earth-Flattening  Approximation  in  the 
Theory  of  Microwave  Propagation,"  Physical  Review,  Vol.  70, 

Nos.  7  and  8,  1  and  15  October  1946,  p.  518. 

5  Pappert,  R.  A.,  and  Goodhart,  C.  L. ,  Waveguide  Calculations  of  Signal 
Levels  in  Tropospheric  Ducting  Environments,  TN  3129,  Naval  Electronics 
Laboratory  Center,  San  Diego,  CA,  25  February  1976. 

6  Budden,  K.  G. ,  The  Wave-Guide  Mode  Theory  of  Wave  Propagation, 

Prentice  Hall,  Ehglewood  Cliffs,  NJ,  1961. 
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such  angles,  which  are  referred  to  as  "eigenangles" .  Whether  the  eigenvalues 
1  are  associated  with  angles  or  with  any  other  physical  parameter,  they  are 

obtained  as  the  roots  of  a  complex  equation  called  a  modal  equation. 
Therefore,  a  major  part  of  any  computation  is  dedicated  to  determining  these 
roots . 


Computations  of  the  lowest  order  waveguide  modes  in  a  duct  environment 
and  the  behavior  of  the  corresponding  fields  were  considered  by  Wait  and 
Spies.  Dresp  developed  a  usable  program  that  computed  field  strengths  in, 
above,  and  below  a  duct.  He  modeled  the  earth  as  cylindrical;  once  a  general 
solution  was  obtained,  he  showed  that  it  can  be  approximated  with  the  aid  of 
the  modified  refractivity  as  the  solution  for  the  flat-earth  case.  Although 
the  formulation  by  Dresp  was  mathematically  elegant,  his  program  was  limited 
to  the  determination  of  20  waveguide  modes.  In  addition,  his  results  were  not 
verified  by  comparison  with  measured  data. 

Pappert  and  Goodhart  (see  Reference  5)  developed  a  code  that  was  verified 

q 

using  measured  data  for  a  surface  duct.  Skillman  and  Woods  found  that  the 
pappert  and  Goodhart  program  predictions  adequately  reflected  measurements  at 
frequencies  below  450  MHz,  taken  in  an  elevated  duct  environment.  However, 
they  found  that  the  Pappert  and  Goodhart  model  failed  to  provide  predictions 
for  comparison  with  measurements  at  2.2017  GHz. 


Wait,  J.  R. ,  and  Spies,  K.  P. ,  "Internal  Guiding  of  Microwaves  by  an 
Elevated  Tropospheric  Layer,"  Radio  Science,  Vol.  No.  4,  April  1969, 
pp.  319-326.  . . ~~ 

8  Dresp,  M.  R. ,  Tropospheric  Duct  Propagation  at  VHF,  UHF,  and  SHF,  MITRE 
Technical  Report  MTR-3114,  Vols.  I  and  II,  MITRE  Corporation, 

Bedford,  MA,  October  1975. 

9 

Skillman,  J.  L. ,  and  Woods,  D.  R. ,  "Experimental  Study  of  Elevated  Ducts," 
Proc.  of  Conference  on  Atmospheric  Refractivity  Effects  Assessment, 
Technical  Document  260,  Naval  Ocean  Systems  Center,  San  Diego,  CA, 

15  June  19  79. 
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Inherent  limitations  can  be  found  in  the  Pappert  and  Goodhart  model.  One 
such  limitation  involves  the  uncertainty  in  any  calculation  that  all 
significant  modes  have  been  found;  that  is,  it  is  possible  in  their  model  to 
miss  eigenvalues.1®  An  attempt  was  made  to  utilize  a  more  effective  search 
method  with  their  model,  but  an  extensive  modification  of  the  mathematical 
formulation  of  the  basic  model  was  required,  making  its  implementation 
impractical . 

Another  limitation  of  the  Pappert  and  Goodhart  model  is  the  requirement 
that  a  "reference  height"  be  specified,  and  the  final  result  is  dependent  in 
some  instances  on  this  reference  height.  Since  the  specification  of  the 
reference  height  requires  some  physical  insight  into  the  particular  problem 
under  consideration,  the  Pappert  and  Goodhart  model  may  be  employed  only  by 
users  with  such  knowledge.  In  this  case,  the  model  is  not  user-oriented.  In 
addition,  no  user's  guide  is  available  for  the  pappert  and  Goodhart  model. 


The  computational  capabilities  noted  above  assumed  a  laterally 
homogeneous  model.  Cho  and  Wait11  performed  work  in  which  only  piecewise 
homogeneity  was  assumed,  and  the  modes  in  each  homogeneous  section  were 
determined.  However,  no  user-oriented  computer  code  for  accomplishing  this 
has  been  developed. 

A  major  obstacle  in  developing  any  computer  program  that  uses  a  waveguide 
model  is  the  calculation  of  all  significant  modes.  The  number  of  modes 
increases  as  the  frequency  and  duct  height  increases.  In  some  cases, 
thousands  of  modes  are  necessary  to  adequately  describe  the  fields;  this 


10  Goodhart,  C.  L. ,  and  Pappert,  R.  A.,  Application  of  a  Root  Finding  Method 
for  Tropospheric  Ducting  Produced  by  Trilinear  Refractivity  Profiles, 
Technical  Report  153,  Naval  Ocean  Systems  Center,  San  Diego,  CA, 

12  September  1977. 

11  Cho,  S.  H. ,  and  Wait  ,  J.  R. ,  Analytical  Study  of  Whispering  Gallery 
Transmission  in  a  Non-Uniform  Tropospheric,  Interim  Report,  Cooperative 
Institute  for  Research  in  Environmental  Sciences,  University  of  Colorado, 
Boulder,  CO,  30  December  1976. 
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further  complicates  the  problem  of  determining  all  significant  modes.  In  the 
hope  of  overcoming  the  problem,  Cho,  et  al.  attempted  to  utilize  a  hybrid 
approach  to  the  problem  in  which  a  waveguide  and  geometrical  optics  series  may 
be  used  to  describe  the  fields.  In  such  a  case,  a  single  geometrical  optics 
term  may  be  used  to  replace  a  vast  number  of  terms  in  the  waveguide  mode 
series.  They  obtained  good  results  for  a  duct  formed  by  an  ideal  refractivity 
profile.  However,  this  method  is  not  applicable  to  more  general  refractivity 
profiles  . 

OBJECTIVE 


The  objective  of  this  work  was  to  develop  a  user-oriented  deterministic 
computer  model  to  compute  the  fields  propagated  in  a  homogeneous  duct 
environment  for  all  frequencies  and  ducts  of  interest. 

APPROACH 

Documentation  was  reviewed  on  existing  computational  capabilities  for 
predicting  propagated  fields  in  a  duct  environment.  Such  capabilities  were 
found  lacking  for  the  frequencies  associated  with  elevated  ducts,  principally 
because  of  the  inability  to  determine  all  the  waveguide  modes  required  to 
describe  the  fields.  After  extensive  analysis,  it  was  found  that  this 
inability  stemmed  from  the  facts  that: 

1.  The  mathematical  functions,  which  are  solutions  to  the  field 
equations  in  the  duct  environment,  became  linearly  dependent  (in  a  numerical 
sense)  for  the  most  significant  values  of  the  argument  of  these  functions;  and 


1  ? 

Cho,  S.  H. ,  Migliora,  C.  G.  and  Felsen,  L.  B. ,  "Hybrid  Ray-Mode  Formulation 
of  Tropospheric  Propagation,"  Proc.  of  Conference  on  Atmospheric 
Refractivity  Effects  Assessment,  Technical  Document  260,  tfeval  Ocean 
Systems  Center,  San  Diego,  CA,  15  June  1979. 
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2.  The  value  of  these  mathematical  functions  could  be  exponentially 
large  or  small,  thereby  exceeding  the  capacity  of  common  high-speed 
computers . 

These  problems  were  overcome  by: 

1.  Expressing  the  solutions  to  the  field  equations  as  linear 
combinations  of  these  functions  in  a  manner  that  assures  linear  independence 
in  a  numerical  sense;  and 

2.  'Utilizing  a  unique  method  for  computationally  expressing 
exponentially  large  or  small  numbers. 

It  was  decided  that  a  Fourier  transform  formulation  of  the  problem  would 
result  in  a  modal  equation  that  was  most  compatible  with  an  available  and 
efficient  eigenvalue  "search"  method.  This  formulation  also  provided  a 
flexibility,  heretofore  unavailable,  for  obtaining  the  field  solution  in  a 
numerically  efficient  manner.  A  computer  code  was  developed  based  on  this 
approach . 

The  adequacy  of  the  mathematical  approach  and  the  corresponding  computer 
model  was  verified  by  comparing  computed  results  with  those  of  other  codes  and 
with  documented  measurements.  Once  this  was  successfully  accomplished,  the 
code  was  used  to  analyze  the  eigenvalues  and  the  fields  for  different  duct 
conf igurations ,  different  frequencies,  and  different  heights  of  the 
transmitting  and  receiving  antennas. 
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ORGANIZATION  OF  REPORT 


The  analytical  and  numerical  basis  of  a  computer  program,  called  DUCT, 
which  has  been  used  sucessfully  in  predicting  propagated  field  strengths  in  a 
DUCT  environment,  is  described  herein.  A  synopsis  of  the  individual  sections 
follows . 

Section  Contents 


2  Mathematical  formulation. 

3  Numerical  procedures  used  to 

determine  the  eigenvalues  of  the 
modal  equation. 

4  Numerical  procedures  utilizing 

the  eigenvalues  to  determine  the 
fields . 

5  Comparison  of  predictions  of  the 

DUCT  program  with  measurements. 

6  Analysis  of  eigenvalue  results 

and  relation  of  specific  eigenvalues 
to  particular  field  contributions. 

7  Analysis  of  results  of  fields  for 

different  parameters  of  interest. 

Appendix  A  Review  of  some  characteristics  of 

the  modified  Hankel  functions  of 
order  one-third. 

Appendix  B  Basis  for  an  efficient  method  for 

evaluating  the  modal  equation,  which 
is  in  the  form  of  a  determinant. 

Section  5  is  independent  of  Sections  2,  3,  and  4  so  that  the  reader  who 


is  interested  only  in  verification  aspects  may  skip  directly  to  Section  5. 
Sections  6  and  7  may  also  stand  alone,  although  Section  3  would  be  beneficial 
for  a  complete  understanding  of  Section  6. 
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SECTION  2 

MATHEMATICAL  FORMULATION 


GENERAL 


In  this  section,  the  problem  of  electromagnetic  propagation  in  an 
atmospheric  duct  is  mathematically  formulated.  This  is  accomplished  for 
horizontally  polarized  radiation  by  using  a  second-order  partial  differential 
equation  to  describe  the  variation  of  a  component  of  the  magnetic  Hertz 
(vector)  potential.  The  solution  of  this  equation  is  obtained  by  utilizing 
Fourier  transform  formalism  to  cast  the  partial  differential  equation  into  the 
form  of  an  ordinary  differential  equation,  and  by  solving  the  latter  equation 
in  terms  of  unknown  coef f icients .  These  coefficients  are  determined  using  the 
boundary  and  radiation  conditions  of  the  system.  The  solution  for  the  Hertz 
potential  is  obtained  in  the  form  of  an  integral  that  is  evaluated  with  the 
aid  of  complex  function  theory  in  terms  of  a  series  of  waveguide  modes. 

The  above  process  is  generalized  to  include  vertically  polarized 
radiation.  Alternative  mathematical  formulations  are  introduced  in  succeeding 
sections  that  will  prove  useful  in  required  numerical  evaluations. 


REFRACTIVITY  PROFILE 


The  index  of  refraction  of  a  medium  is  given  by: 

n(z)  =  /urer  =  /er(z)  (1  > 

where  er  and  ur  are  the  relative  permittivity  and  permeability  of  the 
medium.  It  is  assumed  that  ur  =  1  everywhere.  In  what  follows,  n  is  only  a 
function  of  the  height  z  above  the  ground.  The  refractivity  N  is  defined  as: 

N  =  ( n-1 )  x  1 0^ 


1 1 


(2) 
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and  is  numerically  a  more  convenient  parameter  than  n. 

It  is  often  mathematically  convenient  to  utilize  a  rectangular  coordinate 
system  to  describe  atmospheric  wave  propagation.  In  such  a  formalism,  it  has 
been  found  that  the  earth  curvature  may  be  accounted  for  by  appropriately 
modifying  the  index  of  refraction  of  the  atmosphere.  This  modification  is 
accomplished  by  adding  a  term  to  n  which  would  cause  a  ray  to  bend  in  such  a 
way  that  its  height  above  the  "flat  earth"  at  each  point  would  be  the  same  as 
if  it  were  a  straight  ray  over  a  curved  earth.  The  modified  index  of 
refraction,  m,  is  then  given  by  (see  Reference  2): 

m( z)  =  n( z )  +  z/a  ( 3 ) 

where  a  is  the  radius  of  the  earth  in  the  same  units  as  z. 

The  modified  refractivity  is  defined  as: 

M  =  (m-1 )  x  106  (4) 

The  modified  index  of  refraction  will  have  a  profile  (i.e.,  a  variation 
with  height)  such  as  that  shown  by  the  dashed  line  on  the  right  side  of 
Figure  2.  Tt>  make  the  problem  mathematically  tractable,  this  general  profile 
will  be  approximated  by  a  piecewise  linear  profile  with  L  sections,  such  as 
that  shown  by  the  solid  lines  on  the  right  side  of  Figure  2.  In  the  figure, 

L  =  3  has  been  used.  Each  section  will  represent  an  atmospheric  layer,  or 
region,  the  boundaries  of  which  are  parallel  to  a  flat  earth  located  at  z  =  0. 

The  interface  between  the  ith  and  ith+1  layer  ( 1  _<  i  _<  L-1  )  is  located 
at  z  =  z^,  with  the  layer  i  =  1  closest  to  the  ground.  This  is  illustrated  in 
Figure  2,  for  the  case  L  =  3.  Ehch  layer  is  assumed  to  be  horizontally 
homogeneous,  with  a  modified  index  of  refraction  given  by: 


nr  ( z ) 


z  -H. 

— —  tan  a.  + 
2  l 


1  , 


=  M.(z) 


icf6  + 


<  i  <  L 


(5) 


1  2 
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where  ts  the  value  of  z  at  which  m x  (z)  would  equal  unity,  and  the  slope 
tan  a^/2  is  assumed  small: 

I  tan  a.  I  <<  1  (6) 

l 

From  Equations  5  and  6, 


tn.  (z)  =  1  +  (z-H.  )  tan  a. 

l  li 


1  <  l  <  L 


(7) 


=  1  +  2  M. (z)  x  10 

l 


-6 


The  modified  index  m^(z)  is  assumed  to  be  continuous  across  the  interfaces 
between  the  layers  ( see  Figure  2 ) : 


m^(z^)  =  nu^tz^)  ,  1  <_  i  <_  L-1 


(8) 


MAGNETIC  HERTZ  POTENTIAL 

It  can  be  assumed  that  horizontally  polarized  wave  propagation  is  due  to 
a  radiating  magnetic  dipole  p  =  p  z  oriented  in  the  z -direction  and  located 
at  x  =  y  =  0,  z  =  zT  (see  Figure  2).  In  a  laterally  homogeneous  medium,  the 
wave  propagation  due  to  such  a  dipole  exhibits  axial  symmetry  (about  a 
vertical  axis  through  the  radiating  dipole) ,  and  may  be  obtained  from  the 
z-component  JIz  of  the  magnetic  Hertz  potential  vector  II  .  Thus: 


n  (x,y,z)  =  II  (x,y,z)  z 
z 


(9) 


E  =  -jwy  V  x  ft 
o 


and 


H  =  V  x  V  x  n 
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f  =  the  propagation  frequency, 
UQ  =  the  permeability  in  vacuum, 


and  an  eJ  time  dependence  is  assumed, 


DIFFERENTIAL  EQUATION 


In  each  atmospheric  layer,  II  satisfies  the  partial  differential 


equation: 


V2I[.  +  k  2m.2(z)  n.  =  -p  6(x)6(y)5(z  -  z  ) ,  1  <  i  <  L  (12) 

loi  l  T  —  — 


where 


=  the  value  of  il2  in  the  ith  layer 
=  the  free-space  wave  number  determined  by 
=  the  magnetic  dipole  strength 
=  the  Dirac  delta  function. 


The  Laplacian  operator  v'z  is  defined  as: 


.  2  2 
k  =  <a  p  c 
o  o  o 


V2  =  —  +  —  +  — 
3x2  3y2  3z2 


In  the  ground,  IIZ  satisfies  the  equation: 


V 2 II  +k2n2II  =  0,z<0 
g  o  g  g  - 


(13) 
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where 


rig  is  the  value  of  II  z  in  the  ground  and 

ng  is  the  constant  refractive  index  of  the  ground. 

This  is  given  by: 


where  and  og  are  the  ground  values  for  the  permittivity  and  conductivity, 
respectively. 

A  solution  is  sought  to  Equation  13  in  the  ground  and  to  Equation  12  in 
each  atmospheric  layer,  subject  to  the  boundary  conditions  at  the  ground  and 
at  each  layer  interface  that  the  tangential  components  of  E  and  H  are 
continuous  across  the  boundary.  Using  Equations  7,  10  and  11,  it  may  be 
shown1"*  that  these  conditions  may  be  written: 


"i  +  1 


an. 

i  +  i 
dz 


Z 


z.  1  <  i  <  L-1 

i ,  —  — 


(14) 


(15) 


Tyras , 
Press, 


G. ,  Radiation  and  Propagation  of  Electromagnetic  Waves,  Academic 
New  York,  NY,  1969. 
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An  additional  requirement  on  the  solutions  obtained  is  that  the  radiation 
conditions  be  satisfied.  That  is,  II  must  represent  an  outgoing  wave  as 
z+  -  oo  ,  and  ni  must  represent  an  outgoing  wave  as  z+  +  =°  . 


FOURIER  TRANSFORM  FORMULATION 

The  partial  differential  Equations  12  and  13  may  be  reformulated  as 
ordinary  differential  equations  with  the  aid  of  Fourier  transform  theory  (see 
Reference  13).  The  double  Fourier  transform  of  II  (x,y,z)  is  defined  as: 


n  (u,v,z) 


00  -  j  UX 

=  {  dx  e 

—CO 


oo  -  j  vy 
/  dy  e  II 


(x,y,z) 


(16) 


with  the  inverse  transform  qiven  by: 

1  00  jux  oo  jvy  ~ 

n  (x,y,z)  = - l  du  e  /  dv  e  II  (y,v,z)  (17) 

2 

<  271  ) 

Taking  the  double  Fourier  transform  of  Equations  12  through  15  yields: 


dz 


—  +  k2  [m. 2 ( z ) 

2  o  i 


k2]  f  ni  = 


-p6(z-zT>  ,  1  <_  i  <_  L 


(18) 


2  P 


dz2 

o 

n. 

=  n. 

1 

i 

dll. 

dn. 

i 

i 

dz 

dz 

^  k  2; 


n  =  o 
g 


z  =  z,  1  <  i  <  L-1 

if  —  — 


J 


(19) 


(20) 
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n 

g 


=  n. 


dz 


z  =  0 


where 


2  _ 
P  = 


2 

U 


+  u 


2 


(21  ) 


Thus  solutions  are  sought  for  IL  and  in  Equations  18  and  19,  subject 
to  the  boundary  conditions  in  Equations  20  and  21  and  to  the  radiation 
conditions.  Once  these  solutions  are  found,  they  may  be  used  in  Equation  17 
to  find  the  Hertz  potential  in  each  medium  which,  in  turn,  may  be  used  in 
Equations  10  and  11  to  find  the  electromagnetic  field  vectors. 

SOLUTION  OF  EQUATIONS 

In  Ground 


The  solution  to  Equation  19  is: 


(22) 


where 


y 


k 

o 


2  2 
-<P  Aq  ) 


(23) 


and  Ag  is  a  function  of  the  parameter  p.  n  satisfies  the  radiation  condition 
for  large  negative  values  of  z,  since  it  represents  an  outgoing  wave  in  this 
region.  Tb  assure  that  H  0  as  z  -*■  -“> ,  the  branch  of  y  is  chosen  so  that: 
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Im  Y  <  0. 


(24) 


In  Atmosphere 


The  complete  solution  to  Equation  18  will  be  given  by  the  sum  of  the 
general  solution  of  the  homogeneous  equation  and  a  particular  solution  of  the 
inhomogeneous  equation.  The  homogeneous  form  of  Equation  18  may  be  written: 


+  qi 


ff. 

l 


=  0 


(25) 


where 


q.  ( z)  = 

l 


■<z> 


(26) 


and  Equation  7  was  used.  Equation  25  is  known  as  the  Stokes  Equation  and  its 
solutions  are  given  in  terms  of  Airy  functions  or  in  terms  of  modified  Hankel 
functions  of  order  1/3.  For  purposes  at  hand,  it  is  more  convenient  to 
utilize  the  latter. 

Thus  : 


n. 


i 


A.  h  (q .  )  +  B.h  (q.)  ,  1  <  i  <  L 
ill  l  2  l  —  — 


(27) 


where 


h1  (q) 


(28) 
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(29) 


are  the  modified  Hankel  functions  of  order  1/3  of  the  first  and  second  kind, 
respectively,  and  are  tabulated  in  the  literature.14  hi/3^  and  hi/3^  are 
Hankel  functions  of  order  1/3  of  the  first  and  second  kind,  respectively. 

The  A^  and  in  Equation  27  are  not  functions  of  z,  but  depend  on  the 
parameters  u  and  v.  Their  values  are  determined  from  the  boundary 
conditions.  In  order  to  satisfy  the  radiation  condition  for  large  z  in  the 
Lth  medium: 

A  =  0  (30) 

L 


The  solution  (Equation  27  with  Equation  30)  for  the  homogeneous  form  of 
Equation  18  represents  the  entire  solution  for  every  layer  except  the  one  in 
which  the  transmitting  dipole  is  located.  That  is,  only  for  the  layer 
containing  z  =  zT  is  Equation  18  inhomogeneous.  Assuming  this  to  be  the  Pth 
layer,  a  particular  solution  of  the  inhomogeneous  equation  in  this  layer  is: 

h^ [  qp(z)  )  h2  [  qp  (zT)  ]  ,  z  <  zT 
h2(  qp(z)  i  h1  t  qp  <V  )  ,  Z  >  ZT 


where  W  is  the  constant  Wronskian  for  h^  and  hj  given  by: 


1 4Harvard  Computational  Laboratory,  Tables  of  Modified  Hankel  Functions  of 
Order  One-Third  and  of  Their  Derivatives,  Harvard  University  Press, 
Cambridge,  MA,  1945. 
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w  (hlfh2)  =  tyq)  h2 


(q)  -  h2(q)  hl  (q)  = 


il  f  2- 

IT 


if) 


1/3 


(32) 


and  the  primes  indicate  the  derivative  with  respect  to  the  argument  of  the 
function.  Thus: 


h  (q) 
m 


3h  (q) 
m 

3q 


m  =  1 ,  2 


(33) 


and 


qp'(z) 


% 

dz 


kQ  \  2/3  , 

I  tan  ap|  J 


(34) 


where  Equations  26  and  7  were  used.  That  aquation  31  is  a  solution  may  be 
verified  by  substituting  it  into  Equation  18,  integrating  each  term  from 
zT  -  e  to  zT  +  e  and  letting  z  +  0.  Since  )Tp  is  assumed  finite  at  z  =  zT, 
this  yields: 


dz 


(35) 


which  produces  an  identity  when  Equation  31  is  used. 


Equation  31  may  be  written  more  conveniently  as : 


nF  =  RP  h1  (qP<)  h2  (qP>} 


(36) 
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where 


qP<  =  qp  tmin  (z'zt)] 


qp>  =  qp  Cmax  (z'zt)! 


and 


R 


P 


P _ 

Wq 

P 


It  will  be  noticed  that,  when  P  =  L  and  z  becomes  large,  the  dependence  of 

fl  on  z  is  through  h-[q(z)].  Therefore,  the  solution  of  II  given  in 
P  ^  P/ 

Equation  36,  represents  an  outgoing  wave  as  z  °°,  and  therefore  is 
consistent  with  the  radiation  condition. 


Combining  the  general  solution  (Equation  27)  of  the  homogeneous  equation 
with  a  particular  solution  (Equation  36)  of  the  inhomogeneous  equation  yields: 


II.  =  A.  h  (q.  )  +  B.  h  (q .  )  +R.  h  (q.  )  h  (q.  )  6  ,  1  <  i  <  L  (37) 

l  ill  l  2  l  l  1  i<  2  i>  IP  —  — 


which,  along  with  Equation  30,  is  the  complete  solution  of  Equation  18  in  each 

layer.  £ . „  is  the  Kroneker  delta  function  defined  as: 
iP 

x  _  (  1  -  i  =  P 
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The  coefficients  Ag  in  Equation  22  and  A^  ,  in  Equation  37  must  now 
obtained  from  the  boundary  conditions.  In  the  equations  which  follow,  the 


notation 


q—  is  used  to  indicate  the  value  of  q^  on 


the  boundary  at  z  =  z- 


with  zq  denotinq  ground  level.  Thus: 


q.  .  =  q. (z . ) 

ID  i  D 


Also,  let 


qiT  =  VV 


Substituting  Equations  22  and  37  in  each  of  Equations  21  yields: 


\  =  A1  h1(qio’  +  B1  h 2  +  R1  N  <qi0)  h2(q1T)  61P 


jYA  =  q,  ’  M  ^  (<?0  >  M  {%  >  5  (<*T  > 


from  which  Ag  may  be  eliminated,  resulting  in: 


VV<1,0>  -  6  V’.o"  *  Wqu>>  -  G  h2(q,0>l 

■  -h2lqtT'  V,P[  Vlq,0>  -  Gh,'q,0>1 


where 
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Substituting  Equation  37  in  each  of  Equations  20  yields: 


Vl^ii5  +  Bih2(qii)  +  Rih1(qiT)  h2(qii)  6iP 
=  \+1h1(qi+1.i)  +  Vlh2(qi+1,i) 

+  Ri+lh1(qi+1,i)  h2(qi+l,T)  6i+l,P 


\nr(qii)  +  bi  h2'(qii)  +  Rihi(qiT)  v(qii>  6ip 


V1  \  +  ,  +  Bx  +  1  hVqi+l,i) 

I 

+  Ri+1  h"l(qi+1,i)  h2(qi+1,T;  6i+1,P 


These  equations  may  titf.  recast  into  the  more  convenient  form. 


Aih1(qii)  +  Bih2(qii)  '  Vi  h1(qi+1,i)  '  B  i+1  l‘‘2(qi+1,i) 

=  Ri+1  h1(qi+1ri)  h2(  qi+1,T}  6i+1,P 
-  R. h.  (q.  _)  h  _  (q  .  .  )  6,  „ 
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A.  h' 

i 


!iii 

(qit)  +  si  h2'lqii)  ■  q  { 

qitl' 


+  B  i+1  q/  "2 


*  *1,.  ^ 

-  Rih1(qiT)  h2"(qii)  fiip'  1  <  i-<  L-“  1 


(41b) 


eauations  in  the  2L  unknowns  Ai , 

Equation  41,  39  and  30  nepMsant  2L 
Bi  ,  i  1  1  1  h‘ 

,  the  system  of  equations  which  must  be  solved  to 
As  an  illustration  written  in  the 

obtal„  a*  unsown  co.fticl.nf,  «»•>«  39 

following  form,  assuming  L  =  3  and  P  =  1 • 

A1h1<q10>  +  B1h2<qi0) 

+  B1U2(qH)  ’  A2^(q 


A2h1  (<,21  ) 

-  B2h2lq21) 

-  B2^h2(q21 

A2h1(q22) 

+  B2h2(q22) 

A2h1(q22) 

+  B2h2<q22) 

=  -RTh2(qiT)h1 {q'0] 
=  -R,h1(qiT)h2(qH) 
=  “R1h1(qiT)h2(qH) 


0 

0 


(42) 


where 


<  =  q2/qi  '  q0  =  q3/q2 


and 
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h  (q  )  =  h  '(q  )  -  G  h  (q  )  ,  m  =  1,2  (42a) 

m  1 0  m  10  m  1 0 

If  P  were  not  equal  to  1,  then  the  system  of  equations  would  be  the  same  as 
that  shown  in  Equation  42  except  for  the  right  hand  side  of  the  equations . 
This  will  be  discussed  further  in  Section  4. 

It  is  seen  that  this  system  of  equations  may  be  expressed  in  matrix  forms  as: 

a  £  =  e 

where 

hl (qio' 

|  h1 (q1 1 * 

a  =  |  h/(q11  ) 

0 

0 


(43) 


h2(q10) 

0 

0 

0  1 

h2<qii) 

-  h1  (q21  ) 

-  h2(q21  > 

0 

V(qn> 

-  qR"h1  (q2l  > 

-  qR  h2  (q21  ) 

0 

0 

h1  (q22) 

h1 (q22) 

J 

'  h2(q32>  j 

0 

h," (q22) 

h2  (q22) 

~  qD^  h2*(q32)  j 

r 


£  = 


B,  t 

U  3  J 


and 


B  =  R, 


"  h2(qiT)  h1(q10> 

'  V'V  h2(qi1) 

-  h,(qiT)  h2'(qn: 

0 

0 
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f 


(45) 
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The  solution  for  any  one  of  the  unknown  coefficients  may  be  obtained  through 
the  use  of  determinants.  Thus: 


A. 

l 


Bi 


(47) 


where,  for  1  <_  i  <_  L-1  ,  is  the  matrix  obtained  by  replacing  the  column  of 

a  in  Equation  44  containing  the  coefficients  of  A^  (in  Equation  42)  by  the 
vector  6;  and  for  i  =  L,  |T  |  5  0.  T^,  1  _<  i  _<  L  ,  is  the  matrix  obtained  by 
replacing  the  column  of  a  containing  the  coefficients  of  by  the  vector  0. 

Thus,  for  example: 


"i, 

(q10} 

h2(qi0) 

-  R1h2(qiT)h1(qi0) 

0 

0 

h1 

(qi  i ) 

h2(qi1> 

-  V1(q1T)h2(qi1) 

-  h2(q21} 

0 

TA2  = 

h2<qi1) 

"  R1h1(qiT)h2(qi1) 

’^h2(q21) 

0 

0 

0 

0 

h2(q22) 

-  h2(q32) 

_ 

0 

0 

0 

h2(q22) 

-  qDh2(q32)  _ 

The  notation  |t|  indicates  the  determinant  of  the  matrix  T. 

SOLUTIONS  IN  INTEGRAL  FORM 


Using  Equation  47  in  Equation  38  yields: 


ni(u,v,z)  =  IL  (p  ,z)  = 


lTAi'  W  +  'TBil  W 


+  Rih1(qi<)h2(qi>)6iP  ' 


1  <  i  <  L 


(48) 
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which  may  be  used  in  Equation  17  to  obtain: 

1  «  jyx  »  jvy  ~ 

II  (x,y,z)  =  -  /  dy  e  /  dv  e  H  (p,z)  (49) 

i  2  i 

(  2  TT  ) 

Since  y  and  v  enter  Ih  as  p ,  Equation  49  may  be  written  (References  2 
and  13): 

1  ~  1  ,  "  (2) 
n  (r,z)  =  —  /  p  dp  J  (pr)H  =  —  i  p  dp  H  (irr)  II  (p,z)  (50) 

i  2it0  0  i  4n  C  0  i 


where  C  is  the  contour  in  the  complex  p -plane  shown  as  the  solid  line  in 
Figure  3. 

As  discussed  in  Reference  2,  Equation  50  is  valid  only  when  there  are  no 
singularities  of  IL(p,z)  in  the  first  quadrant  of  the  complex  p-plane,  and 
when : 

/+  P  dp  Hq ^ 1  \pr)  ni  (p  ,z)  -*•  0  ,  (51) 

where  C+  is  the  quarter-circle  contour  of  infinite  radius  lying  in  the  first 
quadrant  of  the  p-plane.  Both  these  conditions  may  be  shown  to  hold. 
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Figure  3.  Contour  for  integral  of  Equation  50 
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SOLUTIONS  IN  TERMS  OF  WAVEGUIDE  MODES 

By  closing  the  contour  C  of  the  integral  in  Equation  50,  the  residue 
theorem  of  complex  variables  may  be  used: 


II.  (r,z)  = 
1 


£  Res  +  -7-  (  /  +  /  )  p  dp  H  *2'(pr)  H.  (p,z)  (52) 

n  n  4tt  v  „  0  l 


where  the  factor  £  Res  is  the  sum  of  the  residues  of  p  H  '  (pr )  JI .  (p  ,  z )  at 
n  n  ox 

the  poles  of  il(p^,z);  the  contour  C~  (see  Figure  3)  is  the  quarter  circle  of 
infinite  radius  lying  in  the  fourth  quadrant;  and  the  contour  B  encloses  the 
branch  cut  which  is  present  as  a  result  of  the  radical  in  the  expression  for  y 
in  Equation  23.  No  other  branch  cuts  are  present  within  the  contour  of 
integration  shown  in  Figure  3  since  the  only  branch  points  of  the  functions 
h1  (q)  and  h2  (q)  are  at  infinity  (see  Reference  14).  The  integral  will 
vanish  over  contour  C~  in  Equation  52  as  in  Equation  51.  The  integral  over 
contour  B  represents  the  surface  wave  contribution  to  the  total  field,  and  is 
assumed  to  be  small  relative  to  the  other  field  contributions.  Therefore,  it 
will  be  ignored  as  well.  Equation  52,  therefore,  becomes: 


II.  (r  ,z)  =  -  i  E  Res 
l  2  n  n 


where  the  residues  correspond  to  waveguide  modes 


An  expression  for  the  residues  in  Equation  53  will  now  be  derived.  Since 
(2  ) 

the  function  H  (pr)  has  no  poles  in  the  complex  p-plane,  the  only  poles  of 

°  (2)  ~  ~ 
the  integrand  p  H  (pr)  n.(p,z)  will  be  those  of  n.(p,z).  But 
o  1  1 

since  h  tq(p)]  and  h  [q(p)]  have  no  poles  in  the  p-plane,  it  is  seen  from 


1 ^Wait,  J.R.,  Electromagnetic  Waves  in  Stratified  Media,  Pergamon  Press, 
New  York,  NY,  1962. 
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Equation  48  that  the  only  poles  of  FL  will  be  those  for  which  |a|  =  0. 
Therefore , 

(2) 

I  O  H. 

(2), 


Res  (p  HA  (pr)  n  (p,z)j  =  Res 
n  0  i  n 


P  Hn(2)(pr)  (IT^I  yy  +  |TB.  |  h2(qj 


or  from  the  well-known  techniques1®  for  evaluating  residues: 

[ITaJ  VV  *  |TbJ  Wjp 


Res  =  p  H  ^ (p  r) 
n  n  0  n 


(54) 


3p 


P  =  P. 


where  p  is  the  nth  value  of  p  for  which  |a|  =  0. 


(2) 


For  all  cases  of  interest,  |pnt|  >>  1,  so  that  Hq  (pnr)  may  be 
approximated  asymptotically  as : 


/ 

(2)  .  .  i tt/4  /  2 

H0  <Pn*>  2  eJ  /^Tr 


-3Pnr 


(55) 


Using  Equations  54  and  55,  Equation  53  may  be  written. 


II.  (r,z) 

i 


1 


Y2TIT 


ej7T/4  Z 
n 


E  e 
n 


-:Pnr 


W  +  ^Bi1  h2(qi) 


-jp  r 
e  n 


P  =  P. 


(56) 


16Churchill,  R.V. ,  Introduction  to  Complex  Variables  and  Applications, 
p.  122,  McGraw  Hill  Book  Co.,  New  York,  NY,  1948. 
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where 


P  =  P 


n 

is  referred  to  as  an  excitation  function,  and: 


Eni  ■  E„i‘vz'V  ■  [lTJ  W  *  !tbJ  w|  <58> 

1  p  =  p 

n 

may  be  referred  to  as  a  height-gain  function.  Hie  functional  dependence  on  z 
occurs  through  h^q^^)  and  h2(q^)  (see  Equation  26).  Hie  parameter  zT  enters 
through  jT^  |  and  |  T Bi  j  . 


ELECTRIC  FIELD  RELATIVE  TO  FREE-SPACE 


Although  the  value  of  E  may  be  obtained  from  Equation  56  by  using 
Equation  10,  it  is  often  more  convenient  to  determine  the  magnitude  of  the 
electric  field  relative  to  free-space: 


where  | E^g |  is  the  magnitude  of  the  field  that  would  be  obtained  at  the  same 
receiving  location  and  using  the  same  source,  but  with  the  propagation  taking 
place  in  empty  space  (i.e.  in  a  vacuum).  Hie  value  of  A  will  now  be  shown  to 
be  proportional  to  JL  for  the  problem  under  consideration. 
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In  free-space,  the  z-component  of  the  magnetic  Hertz  potential  is  given 
by: 


-jV  R  W  q; 

j,  _  P_  e _  _  _P _ P  e _ 

fs  4tt  r  4tt  r 


(60) 


where  Equation  37  was  used  and  q'  is  obtained  from  Equation  34  using  the 
parameters  of  the  actual  environment  (i.e.,  not  the  free-space  environment). 

When  the  Hertz  potential  vector  is  given  by  Equation  9,  Equation  1 0  may 
be  used  to  obtain : 


3IIZ  ^ 
jwu  0 


(61  ) 


in  cylindrical  coordinates.  Using  Equation  60  in  Equation  61  yields: 


J0fs 


k  kip  R  W  q'  ^Qr 
O  P  P  e _ 

4tt  r 


(62) 


where  higher  order  terms  in  1/r  were  neglected.  Applying  Equation  56  to 
Equation  61  results  in: 
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J9i 


—  3  toti 


(2irr) 


'/2 


ej*/4  E 


p  X  E  e 
n  n  n 


'3Pnr 


-jk0oiu 


( 2  it 


r) '2 


ejlt/4  E 


X  E  e 
n  n 


-DPnr 


(63) 


where  it  was  assumed  that: 

p  »  k„  (64) 

n  0 

for  all  p  of  interest.  In  Section  3,  the  approximation  given  in  Equation  64 
n 

is  shown  to  be  valid  for  the  work  described  herein.  By  comparing  Equations  56 
and  63,  it  is  seen  that: 


E0i  =  VWlIi 


(65) 


Defining : 


30  = 


2(27rr)/2  2(2irr) 


]'2 


W  q' 


1/o 

=  2  (  2?rr )  2 


(66) 


and  using  Equations  62  and  65  in  Equation  59  yields: 


A 


4irr  |  II.  | 

1  _  e 

Rp  w  q'  "  P0 


E  A  E  e 


-JPnr, 


n  n  n 


(67a) 
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where  Equation  56  was  used.  A  dipole  strength  of  unity  was  assumed  in 
Equation  66. 

In  Equation  67a,  IL  represents  a  sum  of  complex  numbers  taking  full 
account  of  the  phase  of  each  term.  That  representation  will  be  referred  to  as 
a  "coherent  sum",  a  "mode  sum",  or  a  "vector  sum."  It  is  also  useful  to 
define  the  relative  field  that  would  be  obtained  if  the  phase  of  each  term 
were  completely  random  instead  of  well-defined.  This  will  be  denoted  as  A, 
where : 


A 


=  0 

o 


J E  |X  E  e 
n  n  n 


(67b) 


which  is  similar  to  the  result  that  would  be  obtained  if  the  power 
contribution  of  each  mode  were  added,  rather  than  the  field  contribution. 
Therefore,  the  representation  used  in  Equation  67 a  will  be  referred  to  as  the 
"power  sum"  or  "incoherent  sum." 

In  terms  of  dB  relative  to  free-space,  A  and  A  are  written  as: 


AdB  =  20  1Oq10A 


(67c) 


and 


AdB  =  20  1Ogi0A  <67d) 

In  most  works  on  this  subject  (e.g.,  References  2,  5,  and  8),  the 
function  En  defined  in  Equation  58  is  written  in  the  form: 
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E  =  u { p  ,z)  u(p  ,z  ) 
n  n  n  T 


(68) 


which  demonstrates  the  reciprocity  between  the  transmitter  and  receiver. 

Using  Equation  68  in  Equation  66  would  result  in  precisely  the  formulation 
used  by  Dresp  1 ^  (see  also  Reference  8).  Although  it  is  elegant  in  its 
representation,  the  form  of  Equation  58  is  more  desirable  for  the  numerical 
formulation  of  this  problem. 

VERTICALLY  POLARISED  PROPAGATION 

As  shown  in  Reference  2,  the  magnetic  Hertz  potential,  and  the 

definitions  in  Equations  10  and  11  for  the  electromagnetic  field  vectors,  are 

consistent  with  Maxwell's  equations  when  Equation  1  holds  (i.e.,  when  the 

variation  in  refractive  index  occurs  only  through  a  variation  of  er  but  not  of 

Ur).  For  vertically  polarized  propagation,  on  the  other  hand,  an  electric 
+  (e ) 

Hertz  potential  II  must  be  used,  from  which  the  field  vectors  are  obtained 
through : 


H  =  jeu)  V  x  n 


(e) 


(69) 


and 


-*■  ( e ) 

v  x  v  x  n  '  ' 


(70) 


1 ^Dresp,  M.R.,  Tropospheric  Duct  Propagation  at  VHF,  UHF  and  SHF,  MITRE 
Technical  Report  MTR-3114,  Vols.  I  and  II,  MITRE  Corporation, 

Bedford,  MA,  October  1975. 
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As  shown  in  Reference  2,  these  definitions  are  only  approximately  consistent 
with  Maxwell's  equations.  The  approximation,  however,  is  considered  to  be  a 
good  one  for  the  duct  problems  of  interest. 

-►  (  0  ) 

The  mathematical  formulation  using  the  electric  Hertz  potential  II 
is  identical  to  the  formulation  using  the  magnetic  Hertz  potential,  with  the 
following  exceptions: 


1.  The  magnetic  dipole  strength  p  in  Equation  12  should  be  replaced  by: 
.  (e) 

^  _  -3P 


a)£0n  (zT) 


(  0  ) 

where  p  is  the  strength  of  an  electric  dipole  oriented  with  its  axis  in  the 
z  direction.  This  will  also  affect  the  value  of  Rp  in  Equation  37. 

2.  The  first  of  Equations  15  should  be  changed  to: 

2(e)  2  ( e ) 

n  H  =  n,  (0)  II,  on  z  =  0  (72) 

g  g  i  i 

3.  The  first  of  Equations  21  should  be  changed  to.- 

n  ^  n  =  n  ^(0)  II  on  z  =  0  (73) 

g  g  1  1 

4.  Equation  40  should  be  changed  to: 


37 


ESD-TR-81- 102 


Section  2 


G 


jy 


n  1  2 (0  ) 


2 

n 

g 


(74) 


Of  those  listed  above,  Equation  74  is  the  only  change  that  in  practice 
will  affect  the  calculations.  This  may  be  concluded  from  the  fact  that  8 

o 

(Equation  67)  is  proportional  to  1 /Rp,  and  En  (Equation  58)  is  proportional  to 

Rp  (through  the  dependence  of  |T  |  and  |T^  |  on  R  ).  Since  the  final 

result  is  dependent  on  the  product  B  E  the  value  of  R_  is  therefore  of  no 

on,  f 

consequence.  Also,  Equations  72  and  73  are  manifested  in  Equation  74. 


ALTERNATIVE  REPRESENTATION 


General  Solution  of  Homogeneous  Equation 

The  general  solution  to  Equation  25  may  be  given  in  terms  of  Airy 
functions,  or  in  terms  of  modified  Hankel  functions,  which  are  linear 
combinations  of  Airy  functions.  The  general  solution  to  Equation  25  could 
similarly  be  written  in  terms  of  linear  combinations  of  the  modified  Hankel 
functions.  This  will  be  very  useful  in  later  sections.  Thus: 


n. 

i 


A.  K  (q.  )  +  B. 
Ill  l 


W 


(75) 


would  be  a  valid  general  solution  of  Equation  25,  where: 


w 


C11ih1(qi)  +  C12ih2(qi) 


(76) 


W 


C21ih1(qi)  +  C22ih2(qi) 


(77) 
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The  subscript  1  is  placed  on  the  constants  c  to  indicate  that  there  is  no 

mn 

need  for  these  constants  to  be  identical  in  each  layer.  The  c„_^  would  have 

mm 

to  be  such  as  to  make  the  K^q^  linearly  independent. 

It  is  convenient  to  set: 


C21i  =  °' 


’  22i 


=  1  for  i  =  L 


(78) 


so  that  II  in  Equation  75  would  satisfy  the  radiation  condition  when  Equation 
L* 

30  holds.  It  is  also  convenient  to  define: 


'1  1i 


=  1, 


'  1  2i 


=  0 


'  2 1  i 


=  -e 


4it  j /3 


'22i 


=  1,  i  <  h 


(79) 


and 


c 


1  1  L 


1, 


'1  2L 


e-47ij/3 


(80) 


Substituting  Equations  78  through  80  into  Equations  76  and  77  yields: 


Kl(qi)  =  W 

K2(qi)  =  h2(qi'  ”  e4lr‘i,/3h1  '  1  <  L 


and 

K1(qi)  =  h1(qi)  "  e"4lT:i,/3h2(qi  ) 
K2(qi)  =  h2(qi’ '  i  =  L 


(81  ) 


(82) 
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As  shown  in  APPENDIX  A,  these  definitions  have  the  effect  of  assuring  that  the 
modulus  of  the  two  solutions  K^q^)  and  K2(qi)  are  not  both  exponentially 
large  for  any  value  for  which  Im(q^)  >  0. 

Using  Equations  81  and  82  in  Equation  32,  it  may  be  seen  that  the 
Wronskian  for  K1  and  K2  (as  defined  above)  is  identical  to  the  Wronskian  for 
h1  and  h2: 

w{ki,k2}  =  w|hrh2} 

Particular  Solution  of  Inhomogeneous  Equation 

A  particular  solution  of  Equation  18  was  given  in  Equation  34.  By 
substitution  in  Equation  33,  it  may  be  seen  that  both  of  the  following 
particular  solutions  are  also  valid: 


n 

P 


where  K1  and  K2  are  given  by  Equations  81  and  82.  However,  notice  that 
Equation  85  will  not  satisfy  the  radiation  condition  when  P  =  L. 

Mode  Series  Solution 

Using  the  -representation  rather  than  the  h^-representation,  all  terms 

in  the  mode  series  solution  remain  the  same  as  those  derived  above,  except  h^ 

is  replaced  by  .  thus,  from  Equation  58: 

40 


ESD-TR-8 1-102 


Section  2 


E  = 
n 


W  * 


Bi  1 


W 


P 


(86) 


The  substitution  of  for  would  also  take  place  in  the  expressions 

for  T  . ,  T  .  and  a  in  Equations  44  and  47. 

Al  Bi 
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SECTION  3 

DETERMINATION  OF  THE  EIGENVALUES 

GENERAL 

The  major  computational  effort  required  to  obtain  the  solution  given  in 
Equation  66  is  determined  by  the  eigenvalues  p  =  Pn  which  are  the  zeroes  of 
the  determinant  of  a  (see  Equation  44).  That  is,  the  equation: 


Instead  of  searching  for  solutions  of  Equation  87  in  p-space,  the  search 
may  be  carried  out  for  any  convenient  function  of  p .  In  the  waveguide  mode 
literature,  the  variable: 


is  often  used  (see  Reference  5),  where  0  is  interpreted  as  a  plane  wave  angle 
of  incidence  with  respect  to  the  horizontal  at  a  given  reference  height. 

Dresp  (see  Reference  8)  searches  for  eigenmodes  using  the  variable: 


43 


ESD-TR-81 -1 02 


Section  3 


2 


For  numerical  and  analytical  reasons  which  will  become 
more  convenient  to  search  for  the  roots  of  Equation  87  in  "< 
is,  it  is  more  convenient  to  write  Equation  87  as: 

Gi(qio}  =  |a  Lqion))  I  =  0 

where,  from  Equation  88: 


Thus 


1  - 


l10 


ITT  +  Hi  tan  ai 


0 


tan  a. 


1 


and 


P  =  k. 


1  -H1  tan  a^ 


'1  0 


2/3 


0 


tan  a 


(90) 

clear  later,  it  is 
^-space".  Tfriat 

(91  ) 


(92) 


(93) 


(94) 
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Using  Equation  93  in  Equation  88  yields: 


13 


=  qio  fci  + 


13 


Ui<L,  j  =  i-1,  i 


(95) 


where 


t. 

l 


(|  tan  «1  |  /  |  tan  > 


(96) 


and 


s .  . 
13 


(kQ/|tan  a. 


2^3  ( H  tan  a1 


H.  tan  a.  +  z .  tan  a.  ) 

l  l  3  l ' 


(97) 


The  t^  and  s^^  are  all  real  constants. 


It  should  be  noted  from  Equation  34  that: 


q'  (2) 


(kQ/|tan  a. 


2/3 


tan  a. 
i 


(98) 


is  not  a  function  of  p  (and  therefore  not  a  function  of  q1Q). 

The  roots  of  of  interest  were  shown  to  lie  in  the  fourth  quadrant  of 
the  complex  p-plane.  By  straightforward  conformal  mapping  using  Equation  92, 
the  corresponding  roots  q^^  will  lie  in  the  upper  half  of  the  complex 
q i 0-plane . 
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The  relationship  between  the  and  q1Q  may  be  obtained  from  Equation 

95.  Figure  4  illustrates  this  relationship  for  two  arbitrary  values  of 
q i g  using  a  refractivity  profile  of  the  form  shown  in  Figure  5.  Figure  6 
illustrates  the  relationship  using  the  refractivity  profile  of  the  form  shown 
in  Figure  7. 

NUMERICAL  INSTABILITIES  IN  hi  REPRESENTATION 

It  was  shown  in  the  previous  section  that  the  matrix  elements  in  Equation 
44  may  contain  the  functions  h^  and  h2 /  or  alternatively,  may  contain 
functions  that  are  linear  sums  of  h1  and  h (such  as  the  functions  K1  and  K 2 
in  Equations  76  and  77).  It  can  also  be  shown  that  use  of  the  h1  and  h2 
functions  result  in  numerical  instabilities  of  the  determinant  | ot  |  when 
q10  is  near  the  negative  real  axis  and  has  a  large  magnitude.  Section  6  shows 
that  the  roots  in  this  region  of  the  complex  q^  Q-plane  contribute  most  to 
fields  within  the  duct.  Therefore,  numerical  instabilities  in  this  region 
cannot  be  tolerated. 

From  the  asymptotic  expansions  of  h 1 (q )  and  h2(q)  for  large  |q|,  given  in 
Reference  14  and  APPENDIX  A: 

h2(q)  =  e4lT^3hi  (q)  +  F(q),  <  arg(q)  <  tt  (99) 

For  | q |  large  and  2tt/3  <  arg(q)  <  tt,  h^q)  is  exponentially  large  while  F(q) 
is  exponentially  small.  Thus: 

h2(q)  »  e4lT^3h1  (q)  ,  y-  <  arg(q)  <  ir ,  |q|  >>1  (100) 

If  Equation  100  were  to  hold,  say,  for  q  =  q^  and  q  =  q^  (which  would  be  the 
worst  case  in  branch  B  of  Figure  4),  and  if  Equation  100  were  substituted  into 
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Im  (q) 


Figure  4. 


Relationship  between  .  and  q1Q  in  complex  q-space  for  refractivity 
profile  of  Figure  5.  ^ 
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Equation  44,  then  the  first  two  columns  of  the  matrix  a  would  be  linearly 
dependent.  That  is,  each  term  in  the  second  column  would  be  equal  to  the 
corresponding  term  in  the  first  column  multiplied  by  the  constant 
factor  e41*^3.  This  is  sufficient  to  make  a  singular,  so  that: 

l°t(qio)l  “  °'  a11  qio' 

thereby  making  it  impossible  to  locate  discrete  roots  of  Equation  91.  Of 
course,  the  matrix  a  is  not  singular  in  an  analytical  sense  since  the  term 

F(q) ,  albeit  small,  makes  Equation  100  only  an  approximation.  Nevertheless, 

•• 

the  matrix  a  would  be  singular  in  a  "numerical"  sense,  and  would  preclude  the 

( n) 

use  of  numerical  techniques  to  obtain  the  discrete  roots  q 

Pappert  and  Goodhart  (see  Reference  5)  utilized  a  formulation  employing 
the  functions  h1 (q)  and  h2(q) .  They  avoided  the  numerical  instabilities 
described  above  by  a  "switching"  procedure:  They  argued  that,  for  values  of  q 
for  which  Equation  100  holds,  the  fields  are  evanescent.  Since  the  general 
solution  of  Equation  25  was  shown  to  be  given  by  Equation  27: 


n.  =  A.  h ,  (q.  )  +  B.  h  (q.  ),  1  4  i  <  L,  (101) 

l  ill  i  2  i 


an  evanescent  solution  may  be  obtained  from  Equation  101  by  letting: 
Ai  =  -  Bi.e4lTj/3 


so  that 

Tt.  =  B.  [  h_(q.  )  -  e4lTj/3h  (q.  )|  (102) 

l  l  |  2  l  1  i  1 

which,  from  Equation  99,  may  be  written: 
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n.  =  B.  F(q.  )  (103) 

ill 

As  stated  above,  in  the  region  of  the  q-plane  in  which  Equation  100  holds, 

F(q)  is  exponentially  small  and  therefore  is  an  evanescent  field. 

As  described  by  Fhppert  and  Goodhart  (see  Reference  10),  the  "solution¬ 
switching"  between  Equation  101  and  Equation  102  causes  a  discontinuity  in  the 
function  of  which  the  roots  are  sought  [Gfq^)  in  the  case  at  hand].  In  turn, 
this  could  lead  to  missing  or  duplicating  a  root  already  found. 

CHOICE  OF  ALTERNATIVE  REPRESENTATION 

The  discussion  above  indicates  the  desirability  of  a  solution  in  terms  of 
functions  that 

1.  will  not  be  exponentially  large  simultaneously  in  the  root-search 
region ;  and 

2.  will  reduce  without  "switching"  to  the  form  of  Equation  103  when 
Equation  100  is  valid. 

Functions  that  satisfy  these  conditions  are  the  pair  h^q)  and  F(q)  ,  where 
F(q)  is  defined  from  Equation  99.  Thus: 

n.  =  A.  h.  (q.  )  +  B,  F(q.  )  (104) 

i  i  i  i  i  i 

The  validity  of  using  functions  other  than  h1  and  h^  was  discussed  in 
Section  2. 

In  a  region  where  the  field  is  evanescent,  A^  would  be  expected  to 
approach  0,  thus  leaving  an  equation  in  the  form  of  Equation  103.  That  F(q^) 
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is  exponentially  small  when  |q|  is  large  and  tt  <  arg(q)  <  ti/2  was  noted 
following  Equation  99,  so  that  in  this  subregion  of  the  q10-plane,  h^q)  and 
F(q)  are  not  both  exponentially  large,  h^q)  is  never  exponentially  large  in 
the  remainder  of  the  upper  half  of  the  q1Q-plane  as  shown  in  APPENDIX  A.  Thus 
Equation  104  is  well  suited  for  numerical  determination  of  the  roots  of 
Equation  91.  It  has  the  disadvantage,  however,  of  not  satisfying  the 
radiation  condition.  Since  the  radiation  condition  need  only  be  satisfied  for 
the  solution  in  the  uppermost  region  (i  =  L),  the  form: 

n  =  B  h  (q  )  (105) 

U  u  Z  u  <• 

will  again  be  used  in  this  region. 

The  functions  used  in  the  solutions  given  in  Equations  104  and  105  are 
identical  to  the  functions  K^q^)  and  K2(q^),  as  defined  in  Equations  81  and 
82.  Therefore: 

rtf  =  A.  K1  (q.  )  +  B.  K2(qi)  (106) 

BOUNDS  ON  SEARCH  REGION  -  MODE  ATTENUATION 


As  discussed  previously,  the  eigenvalues  of  Equation  91  will  lie  in  the 
upper  half  of  the  q1Q-plane  (i.e.,  in  the  region  Im(q10)  >  0).  Thus  the  lower 
bound  of  this  region  is  the  real  axis.  Although  a  discussion  of  the  left  and 
right  boundaries  of  the  search  region  will  be  postponed  for  a  later  section 
(Section  6),  an  upper  bound  of  the  search  region  will  be  obtained  here  in 
terms  of  the  maximum  attenuation  of  the  strength  of  a  mode  per  unit  length. 

From  Equation  66,  it  is  seen  that  each  modal  contribution  to  the  total 

-jPnr 

field  has  an  exponential  dependence  e  .  Since  the  p  are  complex  with 

n 

Im(p  )  <  0,  the  amplitude  of  each  mode  will  fall  off  exponentially.  Thus,  if: 
n 

p^  =  a  +  jb,  b  <  0 
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then 

-jp  r  br  -jra 
n  J 

e  =  e  e 

-jPnr 

Since  b  is  negative,  e  will  decrease  in  amplitude  as  r  increases. 

-jpnr 

Typical  values  of  r  are  generally  large  enough  to  make  e  negligible  for 

all  ] b )  greater  than  some  positive  number.  The  relative  decrease  in  the  field 
contribution  from  the  nth  mode  over  a  single  unit  of  length  is  given  by: 


-jp  r( r+1 ) 


-!Pnr 


-JP. 


so  that  the  loss  in  dB  would  be  given  by: 


je  w 

=  -20  log^Q  |e  |  =  -20  log1Q(e  )  =  -20  b  log1Qe  (107) 

Suppose  a  maximum  value  of  is  specified,  say  1^,^.  Then  roots  p^  will  be 
sought,  such  that: 

-b  (20  log,  e)  <  L 

1 0  max 


or 


,  ,  ,  ,  ,  wax 

b  =  Im(p  )  >  —  - - 

n  20  log  e 


(108) 


Equation  108  represents  a  bound  on  the  region  in  p -space  for  which  the  mode 
attenuation  per  unit  length  would  not  exceed  I>max.  This  condition  must  now  be 
converted  to  q1Q-space.  TP  accomplish  this.  Equation  92  is  written: 


q 


10 


(n) 


H  tan  a1  ) 


(109) 
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where 

X 

n 

and 

|X  |  «  1  (110) 

1  n  1 

under  the  assumption  that  |p  |  «  |k^|,  which  will  be  shown  to  be  valid 
later.  Thus: 


where  Equation  110  was  used  to  obtain  the  Thylor  series  expansion  in  Equation 

111. 


Therefore , 


X 

n 


In  particular, 


Im  ( X  ) 
n 


Im  (p  ) 
n 


(112) 


(113) 


Using  Equation  113  in  Equation  109  yields: 


2/3 


i.  -- 


"0 


tar.  a. 


1 


—  lm(p  ) 
*0 


(114) 
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and  using  Equation  1 1 4  in  Equation  108  results  in: 


Im  { q  1 0 )  <  u 


(115) 


where 


U 


2/3 


0 


tan  a 


L 


_ max 

20  logioe 


(116) 


Therefore,  U  as  defined  in  Equation  116  will  be  used  as  an  upper  bound  of 
the  search  region,  with  L  being  specified.  Thus  the  search  for  eigenvalues 

ITlclX 

will  be  carried  out  in  the  band  in  the  q-,Q-plane  defined  by: 


0  <  Im  (qiQ)  <  U 


(117) 


SEARCH  METHOD 


The  procedure  used  to  find  the  zeroes  of : 


Gi(gio)  =  la(qio5 


(118) 


1  8 

in  the  complex  plane  is  described  by  Morfitt  and  Shellman  and  reviewed  in 
Reference  10.  The  procedure  will  be  denoted  herein  as  the  MODESRCH  method. 


1 ®Morfitt,  D.G. ,  and  Shellman,  C.H.,  MODESRCH,  An  Improved  Computer  Program 
for  Obtaining  ELF/VLF/LF  Mode  Constants  in  an  Earth  -  Ionosphere  Waveguide, 
Interim  Report  77T  prepared  for  Defense  Nuclear  Agency,  Naval  Electronics 
Laboratory  Center,  San  Diego,  CA,  1  October  1976. 
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The  method  assumes  that  the  function  G1  is  analytic  everywhere  within  the 
region  of  search,  and  thus  does  not  permit  the  presence  of  poles  in  this 
region.  The  function  G^q^)  satisfies  these  conditions. 

The  MODESRCH  method  was  applied  by  Pappert  and  Goodhart  (see 
Reference  10)  to  the  ducting  problem,  where  their  modal  equation  was  in  the 
form: 

G  (0)  =0  (119) 

where  9  is  related  to  q^Q  and  p  through  Equation  89.  Their  formulation  of  the 
problem  differed  from  the  one  being  presented  here  in  that  the  function  G2  is 
obtained  from  the  fundamental  equation  of  waveguide  propagation: 

G  (0)  =  R  (0 )  R  (0)  -  1  (120) 

2  d  u 

where  R^  and  R^  are  the  reflection  coefficients  "looking  downward"  and 
"looking  upward",  respectively,  from  a  given  reference  height.  Equation  120 
does  not  lend  itself  to  application  of  the  root-finding  method  to  be  discussed 
below  because  it  contains  poles  in  the  region  of  search.  The  function 
G^O),  therefore,  had  to  be  manipulated  to  produce  a  pole-free  function.  This 
fact,  together  with  the  problems  they  encountered  in  solution-switching 
between  the  forms  of  Equations  101  and  102,  detracted  from  the  desirability  of 
applying  this  MODESRCH  method  to  their  formulation. 

Since  the  formulation  of  the  problem  presented  in  this  report  does  not 
suffer  from  the  difficulties  encountered  by  Pappert  and  Goodhart  (i.e.  G^q^) 
is  pole-free  and  no  solution-switching  is  required) ,  the  MODESRCH  method  was 
chosen  to  locate  the  roots  of  Equation  118. 

The  method  utilizes  the  following  fact  from  complex  functional  analysis: 
In  a  finite  region  of  the  complex  q-plane  in  which  no  poles  of  G^ (q)  are 
present  and  in  which  the  only  zeroes  of  G^q)  are  simple  zeroes,  lines  of 
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constant  phase  [i.e.,  arg(G)  =  const]  may  end  only  at  the  boundary  of  the 
region  or  at  a  zero  of  (q) .  Figure  8  (taken  from  Reference  10)  illustrates 
a  rectangular  region  being  searched.  This  figure  shows  that  a  line  AB  of 
constant  phase  on  which  arg(G)  =  0  has  its  ends  at  a  zero  of  the  function  and 
at  the  boundary  of  the  region.  The  line  HG,  on  which  arg(G)=  0,  has  both  its 
ends  on  the  boundary  of  the  region.  It  is  impossible  for  a  line  of  constant 
phase  to  have  each  end  at  a  zero  of  G(q). 

The  MODESRCH  method  thus  starts  from  the  upper  left-hand  corner  of  a 
"search  rectangle"  in  the  q-plane,  and  searches  along  the  left-hand  boundary 
of  the  rectangle  until  a  value  of  q  is  found  at  which  arg[G(q) ]  =  const  =  0  or 
arg(G)  =  180°.  Say  a  phase  contour  arg(G)  =  0  is  found  at  the  point  A,  as  in 

Figure  8.  This  contour  is  then  followed  until  it  either  exits  the  search 

rectangle  or  until  it  ends  at  a  point  within  the  rectangle.  If  it  ends  at  a 
point  within  the  rectangle,  that  point  must  be  a  zero  of  the  function.  Such 
is  the  case  in  Figure  8,  with  the  zero  located  at  point  B.  Since  a  phase 

contour  arg(G)  =  180°  must  also  intersect  that  zero,  this  contour  is  followed 

until  it  exits  the  search  rectanqle  (point  C  in  the  figure).  The  point  at 
which  this  contour  exits  the  rectangle  is  stored,  in  order  to  assure  that  it 
will  not  be  followed  again  later  in  the  search  process.  The  search  then 
resumes  again  at  the  point  A,  and  the  boundary  is  traversed  counter-clockwise 
until  another  value  of  q  is  found  at  which  arg(G)  =  0  or  arg(G)  =  180°.  The 
next  sue)  x)int  in  Figure  8  is  point  C.  But  this  contour  will  not  be  followed 
since  it  has  been  previously  investigated.  Therefore,  the  search  continues 
until  the  contour  arg(G)  =  0  is  reached  at  point  D,  a  zero  is  found  at  point 
E,  and  the  contour  arg(G)  =  180°  exits  the  rectangle  at  point  F.  The  search 
resumes  again  at  D,  and  the  rectangle  is  traversed  until  point  G  is  reached  at 
which  arg(G)  =  0.  This  contour  is  seen  to  exit  the  rectangle  at  point  H 
without  passing  through  a  zero  of  the  function.  The  search  resumes  at  point  G 
and  the  remainder  of  the  rectangle  is  traversed  without  finding  any  additional 
contours  arg(G)  =  0  or  arg(G)  =  180°  which  have  not  been  previously 
investiga ted. 
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Figure  8.  Root  finding  method  for  a  function  G(q). 
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The  procedure  just  described  is  accomplished  numerically  by  dividing  the 
search  rectangle  into  mesh  squares  of  side  length  At  as  shown  in  Figure  8,  and 
by  investigating  the  sign  of  Im(G)  at  the  corners  of  the  mesh.  A  change  of 
sign  between  two  adjacent  corners  indicates  that  a  phase  contour  arg(G)  =  0  or 
a  phase  contour  arg(G)  =  180°  passes  between  these  two  corners.  A  change  in 
sign  of  Re(G)  between  the  corners  of  a  mesh  square  through  which  the  0°  or 
180°  phase  contour  passes,  or  between  the  corners  of  an  adjacent  mesh  square, 
indicates  that  the  phase  contour  arg(G)  =  90°  or  arg(G)  =  270°  is  nearby  and 
thus  a  zero  is  in  the  vicinity.  The  approximate  location  qQ  of  the  zero 
within  a  mesh  square  is  obtained  by  an  interpolation  procedure.  The  precise 

1  g 

location  of  the  zero  is  obtained  using  the  Newton  Raphson  method: 

q.  ,  =  q.  -  Aq.,  i  =  0,1,2,...  (121) 

^1+1  l 


where 


Aq, 


G(qi) 

G^TqTT 


(1  22) 


The  iteration  in  Equation  121  ends  when  the  magnitude  of  the  correction  Aq^ 
to  the  previous  approximation  becomes  less  than  a  specified  small  positive 
number,  i .e .  when: 


I  Aq± i  <  e 


(123) 


1  9 

Pennington,  R.H.,  Introductory  Computer  Methods  and  t&imerical  Analysis, 
Macmillian  Co.,  p.  236  ff.  New  York,  NY. 
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APPLICATION  OF  MODES RCH 


As  mentioned  above,  the  MODESRCH  method  for  determining  the  roots  of  the 
modal  equation  requires  the  definition  of  a  "search  rectangle"  in  the  q1Q- 
plane  where  the  roots  are  sought.  The  upper  and  lower  limits  of  this 
rectangle  were  given  in  Equation  117.  It  is  shown  in  Section  6  that  a  left 
limit  may  be  determined  as  well.  However,  instead  of  defining  the  left  and 
right  sides  of  a  search  rectangle  in  which  all  the  roots  should  lie,  it  is 
convenient  to  perform  the  search  in  subregions.  This  will  have  the  advantage 
of  requiring  less  computer  storage  space  for  the  searching  subroutine.  The 
subregions  are  defined  as  shown  in  Figure  9.  From  Equation  117,  the  lower 
bound  is  at  Im(q1Q)  =  0  and  the  upper  bound  at  Im(q10)  =  U.  [The  MODESRCH 
program  (Reference  18)  automatically  increases  the  search  region  slightly  for 
numerical  reasons.)  The  first  subregion  will  be  the  rectangle  defined  by 
-4  <  Re(qiQ)  <  4.  Thereafter,  subregions  of  4  units  width  are  searched, 
moving  to  the  left  of  the  imaginary  axis,  then  to  the  right,  then  to  the  left, 
etc.,  as  shown  in  Figure  9.  The  search  ends  when  two  successive  regions  are 
encountered  which  contain  no  roots. 


Im(q10) 
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The  size  At  of  the  sides  of  the  mesh  square  in  each  subregion  is  about  .1 

At  *>  .1  (1  24) 


and  the  value  of  e  in  Equation  123  is: 

e  =  At/1000.  (125) 

However,  if,  after  searching  the  first  subregion,  there  is  any  indication  that 
At  is  too  large,  it  is  automatically  made  smaller. 

NUMERICAL  REPRESENTATION  OF  K1  AND  K2  -  EXPONENTIAL  REPRESENTATION  OF  COMPLEX 
NUMBERS 


The  functions  K1 (q)  and  K2(q),  defined  in  Equations  81  and  82,  are 
numerically  evaluated  for  large  values  of  |q|  (i.e.,  |q|  >  4.2)  using  the 

asympotic  expansions  given  in  (Reference  14)  and  summarized  in  APPENDIX  A. 

For  small  values  of  |q|  (i.e.,  |q|  <  4.2)  a  power  series  is  used 

(Reference  14).  In  APPENDIX  A,  it  is  seen  that  for  large  values  of  |q|  the 

magnitudes  of  K.1  and  K2  have  an  exponential  dependence  (see  Equations  A-1 1 
and  A-1 2),  so  that  it  is  possible  for  |kJ  and  |k2|  to  become  exponentially 
large  and  to  exceed  the  upper  numerical  limit  which  a  high-speed  computer  can 
consider.  It  is  this  fact  which  prompted  E&ppert  and  Goodhart  (Reference  3) 
to  "switch"  to  a  solution  of  the  form  of  Equation  102  when  |q|  becomes  very 
large.  Since  the  formula  presented  here  avoids  "switching",  a  different 
scheme  must  be  used  to  permit  evaluation  of  K1  and  K2  even  when  |q|  is 
exponentially  large.  This  scheme  will  be  defined  in  the  following: 

Instead  of  using  two  "words"  to  describe  a  value  of  the  complex  function 
K1  (or  K2)  (i.e.,  a  real  value  and  an  imaginary  value),  three  words  will  be 

used.  These  will  represent  a  real  value,  an  imaginary  value  and  a  real 

exponent.  Thus  for  example: 
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K 


1 


(126) 


where 


K  is  a  complex  number  and 

E  is  a  real  number 

The  value  of  E1  will  be  given  by  one  of  the  exponents  in  Equation  A-1 2.  the 
values  of  j K  )  and  E1  will  be  small  enough  for  a  computer  to  handle  in  a 
single  precision  mode.  This  representation  using  three  numbers  to  describe  a 
complex  number  will  be  referred  to  as  an  "exponential  representation." 

All  arithmetic  using  the  functions  K1  and  K2  will  be  carried  out  without 
evaluating  eE  when  E  is  large.  Thus,  if: 


E 

—  a 
z  =  z  e 
a  a 


-  Eb 
Zbe 


then  the  product : 


z 

c 


z 

a 


z 


b 


will  be  evaluated  by: 


z 

c 


Z*ZH' 

a  d 


E 

c 


E 

a 


+  "b 


(127) 


(1  28) 


(129) 
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4 

The  sum  of  these  numbers  is  obviously  S  =  3e  =  163.79.  However ,  when  summed 
in  the  above  order.  Equation  132  yields: 


S1  = 


+  z  =  (2 


-  -296. 
+  3e  ) 


300 


2e 


300 


(1  34a) 


and 


S  +  S  +  z  =  (2-2)e3°°  =  0  (1 34b) 

1  c 

which  does  not  agree  with  the  actual  result.  The  reason  for  this  disagreement 
is  the  fact  that  a  computer  would  evaluate  2  +  3e  as  2,  thus  losing  the 
second  term. 

The  problem  described  above  may  be  overcome  by  carrying  out  the  sum  in 

En 

the  order  of  decreasing  exponents  and,  if  in  the  partial  sum  =  S^e  the 

value  S  =  0,  setting  E  =  0.  The  example  above  would  then  be  carried  out  as : 
n  n 


S,  =  z  +  z  +  z,  35a) 

1  a  c  b 


so  tha  t 


,  ,  300  „  300  0 

S  =  z  +z  =  (2  -  2)e  -Oe  =0e 
1  a  c 


(1 35b) 


and 
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S 


4  -4 

S  +  z  =  e  (3  +  Oe  ) 


3e 


EVALUATION  OF  THE  MODAL  DETERMINANT 


(135c) 


Determination  of  the  roots  of  Equation  118  required  the  evaluation  of  the 
determinant  of  the  matrix  a  which,  for  L  =  3,  is  given  in  Equation  44  with  h^ 
h2  replaced  by  £2*  Thus: 


Vqio> 

vw 

0 

0 

0 

v«„> 

Vq..> 

-Kt  <q2l’ 

-K2(q2,' 

0 

K2tqi! > 

q2 

-if 

q2  . 

"if  K2lq2>> 

0 

0 

0 

", <q22> 

Vq22> 

~K2(q32) 

0 

0 

K((q22> 

K2(q22» 

q3 

-  -4  K~(q 
q2  2 

(136) 


where 


Km(<W 


K'(qin)  -  G  K  (q  ),  m  =  1,2 
m  1 0  m  1 0 


(1  36a) 


During  the  mode  searching  procedure  the  value  of  |a|  must  be  determined 
numerous  times.  An  efficient  ..iethod  should  be  used  to  accomplish  this. 
However,  the  standard  elimination  methods  require  many  summing  operations 
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which  should  be  avoided  as  much  as  possible  when  expressing  the  K^,  Kj 
functions  in  terms  of  the  exponential  representation  described  above.  The 
evaluation  method  should  also  take  maximum  advantage  of  the  presence  of  zeroes 
in  the  matrix  of  Equation  136. 

APPENDIX  B  describes  a  method  for  accomplishing  this  in  which  the  total 
number  of  multiplication  operations  required  in  the  evaluation  of  the  modal 
determinant  is  on  the  order  of  8(L-1  ),  where  L  is  the  number  of  layers  in  the 
atmosphere.  This  compares  with  (2L-1)!  multiplication  operations  required 
when  evaluating  the  determinant  using  a  cofactor  expansion. 

It  is  also  possible  using  the  method  described  in  APPENDIX  B,  to  avoid 

any  summation  operations  until  all  terms  in  the  sum  have  been  evaluated.  The 

2L—  1 

total  number  of  terms  to  be  summed  will  be  2  ,  and  the  method  described 

following  Equation  134  is  used  to  accomplish  this. 

EVALUATION  OF  THE  DERIVATIVE  OF  THE  MODAL  DETERMINANT 

Part  of  the  MODESRCH  method  for  determining  the  roots  of  G^q^)  in 
Equation  118  involves  the  use  of  the  Newton -Raphson  method  (see  Reference  19) 
described  by  Equations  121  and  122.  Since  Equation  122  includes  the 
derivative  G'(q.|0),  this  derivative  must  be  evaluated. 

If  all  the  elements  a— (q)  of  an  N-by-N  matrix  a  are  functions  of  a 
variable  q,  then  the  derivative  of  the  determinant  of  «  with  respect  to  q  is 
the  sum  of  N  determinants,  each  one  having  the  elements  of  a  different  row 
replaced  by  the  derivative  of  the  elements  of  that  row.  Thus,  for  example, 
for  a  3-by-3  matrix  given  by: 
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the  derivative  of  its  determinant  would  be: 


a11 

ai  2 

3 1  3 

ai1 

31  2 

ai  3 

ai1 

a  1 2 

ai3 

a21 

a22 

a23 

+ 

a21 

a22 

a23 

+ 

a21 

a22 

a23 

3  31 

a32 

a33 

a31 

a32 

a33 

a3; 

*32 

a33 

The  above  rule  for  obtaining  the  derivative  of  a  determinant  may  be 
easily  applied  to  |a|  in  Equation  136.  The  derivatives  of  the  individual 
elements  are  obtained  using  the  formula: 


3  K  (q .  .  )  3  K  (q .  .  )  3  q.  . 

m  l  j  m  i]  13 


3  q 


10 


3 


3  q 


10 


K'  (q.  .  )  t.  ,  1  <  i  <  L, 
mi]  1 


i  -  1 ,  i ,  m 
(1  39) 


where  Equations  95  and  96  were  used  and  the  prime  indicates  differentiation 
with  respect  to  the  argument  q^.  Also: 


3  q 


1  0 


K' (q. . )  =  K  " 
mi]  m 


(q.  .  )  t. 

11  1 


-  q.  .  K  (q.  .  )t. 
l]  m  l]  1 


(140) 


where,  from  Equations  25  and  75,  the  fact  was  used  that  K^fq.^)  is  a  solution 
of  the  equation: 


3  K'(q. . ) 

_ m  1  j 


3  q. 


+  q.  .  K  (q.  .  )  =  0 

l]  m  ni] 


il 


( 1  41  ) 


Equations  139  and  140  may  be  applied  in  Equation  136  to  determine  the 
derivative  of  |a|  with  respect  to  q^  according  to  the  rule  given  in  Equation 
1  38. 
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Each  of  the  individual  determinants  in  Equation  138  will  be  evaluated 
using  the  same  method  as  that  used  to  evaluate  | a  j .  As  previously  discussed, 
each  determinant  is  composed  of  a  sum  of  terms  and,  for  numerical  reasons, 
this  sum  must  be  evaluated  in  order  of  decreasing  exponents  of  the  terms  (see 
Equation  135).  However,  since  the  derivative  of  a  determinant  is  represented 
in  Equation  138  by  the  sum  of  determinants,  it  might  appear  that  the  use  of 
the  summation  scheme  should  be  postponed  until  all  terms  of  all  the 
determinants  in  Equation  138  are  evaluated.  That  is,  in  Equation  138,  each 
determinant  is  composed  of  the  sum  of  six  terms.  lb  avoid  numerical  errors, 
it  wuld  appear  that  the  summation  method  should  not  be  applied  separately  to 
each  set  of  six  terms  representing  the  value  of  each  determinant,  but  rather 
it  should  be  applied  to  the  eighteen  terms  representing  the  entire  summation 
in  Equation  138  (three  determinants  with  six  terms  for  each  determinant). 
Although  this  is  possible  to  accomplish,  it  would  be  extremely  time  consuming, 
particularly  when  the  size  of  the  matrix  becomes  larger. 

The  number  of  terms  in  the  total  sum  representating  the  determinant 
derivative  may  be  made  to  be  the  same  as  the  number  of  terms  in  the 
representation  of  any  individual  determinant  of  which  that  derivative  is 
composed.  (That  is,  in  the  example  of  Equation  138,  the  determinant 
derivative  d|aQ|/dq  may  be  expressed  for  cases  of  interest  as  a  sum  of  only 
six  elements.)  This  stems  from  the  fact  that  the  exponent  in  the  exponential 
representation  of  Km(q^j)  is  equal  to  the  exponent  in  the  exponential 
representation  of  Km'(q^j).  Hi  at  ia,  if: 


_  E  (q. , ) 

K  (q.  .  )  =  K  (q.  .  )  e  m  1]  (142) 

m  lj  mi] 


then 


K' 

m 


(V 


K'(q.  .  )  e 
m  l  j 


E  (q.  .  ) 
m  lj 


(  1  43) 
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with  the  being  the  same  in  both  equations.  It  should  be  noted  that  the 

notation  of  Equation  126  was  followed  so  that  K  (q.  )  is  not  necessarily 

m  l  j  - 

equal  to  the  derivative  of  K  (q  .  .  )  .  The  fact  that  the  E„  are  the  same  in 
H  m  13  m 

Equations  142  and  143  may  be  seen  from  examining  the  derivatives  of  the 
asymptotic  expansion  of  Km  in  APPENDIX  A. 

Thus  in  each  of  the  determinants  of  which  9|o  |/3q1  is  composed,  the 
value  of  the  exponent  of  the  (i,j)  term  would  be  identical.  Hiis  fact  will 
lead  to  the  terms  of  the  sums  representing  each  determinant  having  the  same 
exponents  for  each  determinant.  Tb  illustrate  this.  Equation  138  will  be  used 
wi  th: 


E.  . 

a  =  ~  e 

ij  ij 


(144) 


(1  45) 


Then,  from  Equation  138: 

d,aol 


.  =  S ,+  +  S, 

dq  12  3 


where 


G1  C2  C3  £4  £5  £6 

S1=b1,e  +b126  +b13S  *b14S  +bl5S  +  b16e 


£1  £2  £3  £4  £5  £6 

S2  =  b2l@  +  b226  +  b236  +  b24S  +  b25S  +  b26e 


S 3  =  b3ie  1  +  b326  2  +  b336  3  +  b346  '  +  b35S  '  +  b36S  &  <146) 
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where  the  e^'s  are  the  same  for  and  S^,  1  <  i  <  6. 

Since  these  s^'s  are  the  same,  the  respective  terms  may  be  summed  without 
incurring  numerical  difficulties.  Thus: 

d|“0l  S  £2  £3  £4  £5  £6 

=Cle  +  c2e  +  c3e  +  c^  +  c5e  +  c  (147) 


where 


=  b1i  +  b2i  +  b  3i  ' 


1  <  i  <  6 


(  1  48) 


From  Equations  146  and  147  it  is  seen  that,  using  the  summation  scheme  for 
exponential  representation,  the  number  of  terms  to  be  summed  in  order  to 
obtain  the  determinant  derivative  is  the  same  as  the  number  which  must  be 
summed  for  the  determinant  itself. 


SAMPLE  RESULTS 


A  method  has  been  described  above  for  numerically  determining  the  roots 


L10 

refractivity  profile  shown  in  Figure  10,  the  results  for  the  case  of 
transmission  frequency  of  2.2017  GHz  are  shown  in  Figure  11.  The  value  of 
Lmax  used  in  Equation  107  was  .375  dB/km. 


Each  "dot"  in  the  figure  represents  a  mode.  For  comparison  purposes,  an 
illustration  of  the  same  modes  in  a  space  similar  to  one  utilized  in 
Reference  7  is  also  included  as  Figure  12.  Figure  12  plots  the  location  of 
the  roots  in  terms  of  Ret  9)  with  9  given  in  Equation  89,  and  the  mode 
attenuation  per  kilometer  defined  by  Equation  107. 
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In  Figures  1 1  and  1 2,  the  corresponding  regions  in  which  the  roots  lie 
are  indicated.  Thus  the  roots  lying  in  region  A  of  the  q^-plane  in  Figure  11 
correspond  to  the  roots  lying  in  region  A  of  the  complex  plane  shown  in 
Figure  12.  The  approximate  number  of  roots  in  each  region  are: 


Region 

A 

B 

C 

D 


No.  of  roots 

34 

35 
67 

344 


An  interesting  result  in  the  comparison  of  the  two  figures  is  the  large 
variation  in  root-spacing  that  occurs  from  region  to  region  in  Figure  12. 

Thus  the  34  roots  in  region  A  in  Figure  11  are  essentially  located  at  a 
"point"  on  the  scale  of  Figure  12.  Similarly,  the  67  roots  in  region  C  in 
Figure  11  are  all  clustered  in  a  relatively  small  region  of  Figure  12.  The 
magnified  plot  of  the  roots  shown  in  Figure  13  (corresponding  to  Figure  12), 
demonstrates  this  to  be  true.  The  reason  for  this  behavior  is  that  the 
abscissa,  Re (6),  is  an  entirely  different  quantity  from  the  ordinate, 
attenuation/km .  Indeed,  had  the  roots  been  plotted  in  a  6-plane  in  which  the 
axes  were  Re ( 9 )  and  Im  (9),  then  the  behavior  would  have  been  appreciably 
better.  Nevertheless,  it  is  to  be  noted  that  this  poor  behavior  was  not 
observed  in  similar  plots  presented  in  Reference  9.  this  is  probably  due  to 
the  difference  in  mathematical  formulation  in  Reference  9,  including 
specifications  of  a  "reference  height"  within  the  duct  (see  Section  5). 


It  is  apparent  from  Equation  89  that,  for  all  the  roots  of  interest,  the 
corresponding  values  of  are  very  close  to  the  value  of  kQ,  because  of  the 
smallness  of  the  values  of  9  for  the  roots  in  Figure  12. 
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Eigenvalues  obtained  in  0-space  using  refractivity  profile  of  Figure  10. 
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NORMALIZATION  OF  REFRACTIVITY  PROFILE 


It  will  now  be  shown  that  the  roots  q1Q^n'  obtained  for  a  given 
refractivity  profile  M(z)  are  identical  to  the  roots  that  would  be  obtained 
for  a  profile  with  the  same  refractivity  gradients,  but  normalized  so  that 
M(  0 1  =0.  That  is,  the  roots  would  be  the  same  for  the  two  refractivity 
profiles  shown  in  Figure  14. 

Using  Equation  7,  let  the  unnormalized  modified  refractivity  be  given  by 


2  M. (z)  .  10“6  =  (z  -  H. )  tan  a. ,  1  <  i  <  L  (149) 

1  li 

Then  a  modified  refractivity  profile  M^(z)  with  the  same  gradients 
(characterized  by  tan  )  would  be  given  by: 


2  M  (z) 
1 


10 


-6 


(z  -  H  )  tan  a. ,  1  <  1  <  L 

1  1 


(1  50) 


Since  it  is  desired  that  M  (z)  have  the  value  0  when  z  =  0: 

l 

H  i  0  (151) 

l 

Since  M^(z)  and  M  (z)  have  ti.e  same  gradients,  the  difference  A  between 
—6  -  _  ^ 

2M  (z)  x  10  and  2M  (z)  x  10  at  each  value  of  z  will  be  constant.  Thus: 
i  i 

__  IT  tan 

M.  (Z)  -  M~(z)  =  const  =  M.  (0 )  -  M.  ( 0 )  =  M. (0)  — - 1 - 

1  1  1  1  1  2x1  O'6 

(152) 


_ A _ 

2x1  0 
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Figure  14.  Illustration  of  an  unnormalized  and  a  corresponding  normalized 
refractivity  profile  . 
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or,  using  Equations  149  and  150  in  Equation  152, 

(Z  -  H. )  tan  a.  -  (Z  -  H. )  tan  a.  =  -  H  tan  a. 
ii  ii  11 

or 


H,  tan  a.  -  H.  tan  a.  =  -  H.  tan  a.  (153) 

1  1  i  i  i  i 


Now  the  eigenvalues  of  the  modal  determinant  for  the  refractivity  profile 
M^(z)  will  be  the  same  as  the  eigenvalues  of  the  modal  determinant  for  the 
refractivity  profile  M^(z)  if  the  q  _  ,  1  _<  i  £  L,  j  =  i-1,i  are  the  same  for 
each  profile.  For  the  profile  M^(z)  given  by  Equation  149,  the  values  of  q^j 
will  be  given  by  Equations  95  to  97,  with  Equation  97  representing  the 
dependence  of  q^  on  the  as: 


s  .  . 

13 


tan  a  -  H.  tan  a  +  z.  tan  a. ] 
1  i  13  1 


(154) 


For  the  profile  Nh(z)  given  by  Equations  150  and  151,  the  values  of  q^j  will 
be  given  by  Equation  95  to  97  but  with: 


s  .  . 
13 


2/3 


0 


tan  a. 


(-  H.  tan  a,  +  z.  tan  a.] 
1  13  1 


(155) 


But  from  Equation  153,  the  value  of  s^j  in  Equation  155  is  the  same  as  the 
value  of  s^j  in  Equation  154.  Therefore,  the  values  of  q—  are  the  same  for 
the  two  profiles,  and  thus  the  values  of  q^Q  are  also  the  same. 

It  is  emphasized  here  that  the  fact  that  the  values  of  q^^0^  are  the 
same  for  the  normalized  and  unnormalized  profiles  does  not  imply  that  the 
fields  obtained  in  Equations  67  for  the  two  profiles  would  be  the  same.  In 
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general,  these  fields  would  not  be  the  same  for  the  two  profiles  since, 
although  the  q1Q(n)  are  the  same,  the  pn  are  not.  The  pn  affect  the  value 

“jpnr 

of  e  and  A  in  Equations  67. 
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SECTION  4 

NUMERICAL  DETERMINATION  OF  THE  FIELDS 
RELATIVE  FIELDS  IN  TERMS  OF  qiC)(n) 

Once  the  zeroes  of  the  modal  determinant  of  the  problem  are  found,  the 
electric  field  magnitude  relative  to  free  space  due  to  a  transmitter  at  height 
zT  may  be  determined  from  Equation  67a  at  a  height  zR  and  distance  r  from  the 
tr  an  s.;ai  tter .  Thus 


A 


l  X 


n  n 


-jP 

E  e 
n 


r 


n 


(1  56) 


where,  from  Equations  57,  58  and  66: 


6q  =  2  /  2irr 


(157) 


X  (p  ) 
n  n 


(158) 


E 

n 


(pn' 


V 


^TAR(qT} 


w  + 


|TBR(qT)ik2(V 


Mo 


( 1  59) 


M  0 


(n) 


In  Equation  156,  pn  is  obtained  from  the  eigenvalue  q1(^n'  through 
Equation  94: 
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/ 

P  =  k  J  1  -  H  tan  a 


q 


( n) 

1  0 


I  tan  aj 


2/3 


(160) 


Equation  158  may  be  expressed  in  terms  of  q 


10 


(n) 


by  writing: 


AT 


■  /dial 

3  Ml 

[V3qio 

3p  j 

P  =  Pn'  qi0  =  qi0 


(n) 


But,  from  Equation  92, 


3q 


1  0 


3p 


2p 


2/3 


'0 


tan  a 


Therefore 


X  = 
n 


2/3 


(1  61  ) 


2  / p 


0 


tan  a 


9  |ct  I 
3qic 


qio  =  qio 


(n) 


The  value  of  in  Equation  161  may  be  obtained  from  Equation  160  and  the 
value  of  3 |  ot  | /3q 1 ^  may  be  obtained  using  the  method  discussed  by  Equations 
137  to  148. 


The  factors  X^,  E^  and  exp  (-jp^r)  in  Equation  156  are  evaluated  using 
the  exponential  representation  illustrated  in  Equation  126,  and  the  product  of 
these  factors  is  evaluated  using  Equations  127  to  129. 

Referring  to  Equation  159,  the  qR  are  defined  using  Equations  95  to  97 
by: 


qR  = 


qi0tR 


+  S. 
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where 


and 

2/3 

|hi  tan  a1  -  Hr  tan  aR  +  zr  tan  aR  J  (164) 

where  zR  is  the  height  of  the  receiver  located  in  region  R.  If  the 

transmitter  is  at  a  height  zT  which  is  located  in  region  T,  then  q,^  is  given 

by  Equations  162  to  164  with  R  replaced  by  T.  The  values  of  qR  and  qT  used  in 

Equation  159  are  those  for  which  q10  =  q1Q(n)  in  Equation  162. 

DEFINITIONS  OF  THE  MATRICES  AND  TRR 

Equation  159  requires  the  values  of  the  determinants  | T  |  and  | T  |  . 

Aa\.  dR 

The  matrices  T  and  T  will  first  be  considered.  These  were  discussed 

AK  BK 

briefly  following  Equation  47.  For  clarity,  they  will  be  redefined  here  using 
the  solutions  rather  than  the  solutions. 
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Tlie  modal  matrix  for  a  three-layer  atmosphere  is  given  by 
(see  Equation  136): 


A/So1 

Vqio> 

0 

0 

'  Vq„> 

K2(qt 1 1 

-Vq21) 

-K2(q2,( 

a  = 

q2 

"if  Vq21> 

'J2 

-  5^  K2<q21 > 

\  0 

0 

K1(q22) 

K2(q22> 

\o 

0 

K1(q22) 

K2<q22> 

■K2(q32) 


qT  K2(q32)- 
(165) 


where,  as  in  Equation  136a, 


K  (q  )  =  K'(q  )  -  G  K  <q  )  ,  m  =  1,2 
m  1 0  m  1 0  m  1 0 


(166) 


As  in  Equation  46,  a  vector  BT  is  defined  which  represents  the  contribution  to 
the  boundary  conditions  of  the  particular  solution  of  the  inhomogeneous 
differential  Equation  18.  The  subscript  T  indicates  the  region  in  which  the 
transmitter  is  located  (1  <  T  <  L) .  For  the  example  of  L  =  3, 


(1  67) 


(1  68) 
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The  values  of  the  8^  will  be  given  below.  Thus,  when  the  transmitter  is  in 
region  1,  Equation  167  is  used;  when  the  transmitter  is  in  region  2,  Equation 
168  is  used;  and  when  the  transmitter  is  in  region  3,  Equation  169  is  used. 

The  matrix  TftR  is  constructed  by  replacing  the  first  column  of  the  Rth 
pair  of  columns  of  a  by  the  vector  0T,  1  R  <_  L-1,  and  the  matrix  TBR  is 
constructed  by  replacing  the  second  column  of  the  Rth  pair  of  columns  of  a  by 
0T,  1  £R  <_L-1.  The  matrix  TftL  is  not  considered  as  explained  following 

Equation  47.  The  matrix  ?BL  is  constructed  by  replacing  the  last  column  of 
a  by  the  vector  BT. 

Thus,  TA1  is  formed  by  replacing  the  first  column  of  a  by  6T;  and 
Tft2  is  formed  by  replacing  the  third  column  of  a  by  BT>  TR1  is  formed  by 
replacing  the  second  column  of  a  by  0T;  TB2  is  formed  by  replacing  the  fourth 
column  of  a  by  8T;  and  TR3  is  formed  by  replacing  the  fifth  column  of  a  by  0T* 

A  discussion  of  the  elements  of  the  0T  vector  follows.  The  elements  of 
0T  for  T  =  1  were  given  in  Equation  46  using  the  solutions  in  terms  of  the  h^ . 
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Suppose  the  solution  of  the  inhomogeneous  differential  equation  in  region 
1  were  given  in  terms  of  Equation  84: 


S,  -  R,  K,(q,<)  K2(q,>)  (.70, 

where  q1<;  and  q1  >  are  defined  following  Equation  36.  Then  the  elements  of  B1 
may  be  given  in  terms  of  the  solutions  as: 


0 

6 

0 


1 1 

12 

13 


R,  IC2"»,Tll',<q,0) 


R1  K,  "t1T)K2<'I1  1  ’ 


(1  7.  > 


where  the  01m  refer  to  the  corresponding  elements  in  the  vector  of  Equations 
167.  If  the  solution  were  given  by  Equation  85: 


n,  -  -  R,  Vq^Wjlq,,) 

then  the  following  set  of  values  of  the  0lm  would  be  valid: 


611 

=  R1 

VhT’Vho’ 

S1  2 

=  R1 

K2<qiTlK,<qi, 1 

(173) 

S13 

-  R1 

K2lq,TIK1lq,,» 

In  an  entirely  similar  manner,  if  the  solution  in  region  2  is  given  by: 
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"2  '  *2  WW 


then  the  elements  of  82  are: 


62,  '  E2  WV^,’ 


>22  ‘  *2  5: 


623  ‘  -  “2  ^1^21^2^22^ 


624  ‘  -  R2  R1  ^2T ^2*^22 ’ 


and  if  the  solution  is  given  by: 


n2  '  -*2  V^’VW 


then 


S2,  -  -  R2  *S  ^2T  ^R2  '^21  * 


*22  ’  -  R2  5^  > 


B23  ’  R2 

S24  *  R2  K2<q2T)KT<<:,22) 


In  region  3,  if  the  solution  is  given  by: 


n3  =  R3  K,(<?3<)K2(q3>) 


(174) 


(175) 


(1  76) 


(1  77) 
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then  the  elements  of  83  are: 


<*31  *  *3  K2<q3T>K,(C132l 
®32  R3  K2(q3T,Kl (q32! 


and  if  the  solution  is  given  by: 


"3  '  -  S3  Vq3>,K2!q3<’ 


the  elements  of  0^  are: 


(179) 


(180) 


31 


R3  K1<q35)K2(q32> 


S32  '  R3  q'  K1(q35)K2(q32) 


(181) 


In  region  3,  the  solution  II ^  must  satisfy  the  radiation  condition,  thus 
precluding  the  use  of  Equations  180  and  181.  In  region  3,  therefore,  83  will 
be  defined  using  Equation  179. 


Theoretically,  either  Equations  171  or  173  may  be  used  to  define 
8^;  and  either  Equations  175  or  177  may  be  used  to  define  8  Numerically, 
however,  there  is  often  a  strong  preference  as  to  which  expression  for  8-j  and 
B-,  to  choose.  The  basis  of  the  choice  will  be  somewhat  similar  to  the  basis 
used  to  choose  Equation  179  instead  of  Equation  181  to  define  the  vector  83  in 
region  3.  Just  as  in  region  3  the  form  of  the  solution  is  used  which  "best 
satisfies"  the  boundary  and  radiation  conditions,  so  too  the  form  of  the 
solution  which  most  closely  characterizes  the  field  in  each  respective  region 
will  be  used.  Heuristically ,  the  general  solution  will  have  to  "work  harder" 
to  satisfy  the  boundary  conditions  when  the  particular  solution  is  farther 
from  the  exact  solution.  Conversely,  the  general  solution  will  be  more 
numerically  correct  if  the  particular  solution  more  closely  characterizes  the 
field  in  the  region  of  interest. 
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lb  determine  which  solution  best  characterizes  the  field  in  a  given 
region,  the  following  general  information  will  be  used,  and  verified  in 
Section  7.  In  a  duct  environment, 

1.  The  only  instance  in  which  the  field  in  region  T  could  be  much  less 
at  the  upper  boundary  of  the  region  than  at  the  height  at  which  the 
transmitter  is  located,  is  when  tan  aT  <  0,  T  >  1  and  the  field  will  be  that 
due  to  the  "trapped  modes"  (see  Section  6). 

2.  "Trapped  modes  are  those  for  which  the  eigenvalues  lie  near 

the  negative  real  axis,  which  implies  from  Equations  162  through  164  that 
qT<^n^  and  qT>^n'  near  the  real  axis. 

From  Equations  162  through  164: 

Re(qT<)  <  Re(qT>)  when  tan  aT  >  0  (182) 


Re(qT<)  >  Re(qT>)  when  tan  ctT  <  0  (183) 

When  Re(q.|Q)  <<  0,  it  may  be  seen  from  APPENDIX  A  and  Equations  182  and  183 
that : 

I K-1  (qT< )  I  >  lKi(qT>^'  tan  aT  >  °  When  lKi(qT<)l  >>  1  (184) 


I K !  ( qT> )  I  >  liS(qT<)l'  tan  aT  <  0  when  lKi(qT>)l  >>  1  (185) 
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when 


lK2(qT<) 


'K2(qT>) 


(1  86) 


From  Equations  184  to  186,  for  tan  aT  >  0,  it  may  be  concluded  that  the  field 
in  region  T  due  to  trapped  modes  will  not  be  appreciably  less  at  the  upper 
boundary  of  the  region  than  at  any  point  within  the  region  when  the  solution 
in  that  region  is  given  by: 


II  ~  K  (q  )K  (q  ),  tan  a  >  0 
T  1  T>  1  T<  T 


(1  87) 


For  tan  aT  <  0,  the  field  may  be  appreciably  less  at  the  upper  boundary  of 
region  T  than  at  any  other  point  within  the  region  when  the  solution  in  that 
region  is  given  by: 


nT  ~  V'WV'W'  tan  aT  <  0  (188) 


Therefore,  from  item  (1),  the  particular  solution  of  the  differential 
equation  will  be  given  by: 


f  Ki(qT<)K2(qT>)'  tan  aT  >  0  or  T  =  1 


( 1  89a) 


=  R„ 


< 


~  Ki (qT>)K2(qT<)'  tan  aT  <  0  and  T  *  ' 


(1  89b) 


88 


ESD-TR-81 -1 02 


Section  4 


The  values  of  the  elements  in  the  vector  8T  will  be  determined  based  on 
whether  the  solution  is  given  by  Equation  189a  or  by  Equation  189b.  Thus,  the 
elements  of  6(  will  always  be  given  by  Equation  171.  The  elements  of  gj  will 
be  given  by  Equation  175  if  tan  >  0  and  by  Equation  177  when  tan  a <  0. 

All  evaluations  of  the  | T  ]  and  | T  | ,  as  well  as  of  E  in  Equation  159, 

AR  dK  “ 

are  carried  out  using  exponential  representation. 


EVALUATION  OF  THE  DETERMINANTS 


V  MD  'V 


The  form  of  the  matrices  T  and  T  (i.e.,  the  locations  of  their  zero 

AR  BR 

elements)  depends  on  the  value  of  R  and  the  value  of  T.  The  value  of  T  would 

affect  whether  8  is  given  by  Equation  167,  168  or  169.  The  value  of  R  would 

affect  the  column  of  a  (Equation  165)  which  the  vector  6  replaces  to  form 

T  or  T 
AR  BR . 


It  is  clear  from  the  three-layer  example  used  in  Equations  165,  167,  168 
and  169  that,  if  R  =  T,  the  zero  elements  of  T  and  T  are  in  the  same 

A  R  BR 

location  as  the  zero  elements  of  a.  The  same  algorithm  used  to  evaluate  |a| 
(see  APPENDIX  B)  will  then  also  be  applicable  to  the  evaluation  of 


AR' 


and 


BR' 


Tb  demonstrate  the  general  evaluation  procedure  when  R  $  T,  only  a  single 
example  will  be  used,  with  the  application  to  other  cases  being 
straightforward.  For  this  example,  T  =  1,  R  =  2  will  be  used,  and  the 
determinant  |Ta2|  will  be  evaluated.  In  this  case,  the  matrix  Tft2  will  have 
the  form: 
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o>  Vqio> 


Vqn>  W 


1  1 


1  2 


0 


\ 


\ 

o  \ 


A2 


K'(qu) 


13 

0 


-  ST  K2("2, ' 


V>22> 


K2(q22> 


K2(q32> 


K2«q32> 

(190) 


Now  comparing  T  with  a,  it  is  seen  that  has  non-zero  elements  in  the 

same  locations  as  the  non-zero  elements  of  a,  except  T  has  the  element 

at  a  position  in  which  the  matrix  a  has  a  zero.  But  from  the  definition 
of  a  determinant  in  terms  of  a  cofactor  expansion,  it  may  be  shown  that: 


ita2i  =  lTA2| 


+  T 


A2 


(191  ) 


1  1 


B12  =  613  =  ° 


Now  T 


A2 


has  non-zero  elements  in  the  same  locations  as  the  non-zero 


611=  0 


elements  of  |a|,  and  may  therefore  be  evaluated  using  the  same  procedures  as 
those  used  to  evaluate  |a|  (APPENDIX  B) .  The  second  term  on  the  right  side  of 
Equation  191  may  be  written  as  the  product  of  0^  and  its  cofactor: 


A2 


12 


=  6 


B,,l«(B,r 


(192) 


1  3 


where  M(  01  1  )  indicates  the  matrix  which  is  the  minor  of  0^.  and  is  given  by: 


90 


r 


ESD-TR-81 -1  02 


M(  8 


1  1 


1  02 

Section  4 

/ 

\ 

/K1(qi1> 

K2<qu’ 

-  K2<q2,' 

0  \ 

1 

K?q,  1> 

q2 

-  51  VO 

0 

\  ° 

\ 

0 

VO 

-  / 

\  0 

0 

k;  (q  ' 

q3  „  / 

-  -4  K'ICq  )  / 

(193  ) 

2  ^22 

q2  2  32  ' 

The  determinant  of  M(811)  (also  called  the  cofactor  of  B^)  is  just  the 
product  of  the  determinants  of  the  -upper  left  2-by-2  matrix  and  the  lower 
right  2-by-2  matrix: 


VO 

VO 

M(3n)| 

|  = 

X 

vo 

Vqn> 

K2(q22) 


K2(q22) 


K2(q32) 


32 


(190 


The  evaluation  of  the  determinants  on  the  right  side  of  Equation  194  is 
straightforward.  However,  the  first  determinant  may  be  identified  as  tne 
Wronskian  given  in  Equation  83,  so  that  this  determinant  need  not  be 
numerically  calculated. 

Substituting  Equation  194  into  Equation  192,  the  result  may  be  used  in 
Equation  191  to  evaluate  |T  |  ’ 
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SECTION  5 

VERIFICATION  OF  CALCULATION  METHOD 


GENERAL 

A  method  is  presented  in  the  previous  sections  for  calculating  the  fields 
relative  to  free  space  in  a  duct  environment.  The  method  has  been  programmed 
in  a  computer  model  called  DUCT.  In  this  section,  the  predictions  using  DUCT 
are  compared  with  measured  fields  in  an  environment  containing  surface  ducts 
and  in  an  environment  containing  elevated  ducts.  In  each  case,  the  measured 
fields  are  compared  with: 

(1 )  the  mode  sum  predictions  (Equation  67c)  in  which  the  phase  of  each 
modal  contribution  is  fully  taken  into  account;  and 

(2)  the  power  sum  predictions  (Equation  67d)  in  which  the  phase  of  each 
modal  contribution  is  assumed  to  be  random. 

SURFACE  DUCTS 


Measurements  in  an  environment  containing  a  surface  duct  were  documented 
in  several  references  '  (see  also  Reference  5).  The  refractivity  profile 
used  to  calculate  the  fields  is  illustrated  in  Figure  15.  It  is  the  same  as 
that  used  by  Pappert  and  Goodhart  as  an  approximation  to  the  profile  observed 
over  the  propagation  path  during  the  period  in  which  the  measurements  were 
made.  It  is  normalized  to  zero  at  the  ground.  In  Section  3,  it  was  shown 
that  this  normalization  does  not  affect  the  location  of  the  eigenvalues 


20 

Pappert,  R.A.,  and  Goodhart,  C.L.,  "Case  Studies  of  Beyond-the-Honzon 
Propagation  in  Tropospheric  Ducting  Environments,"  Radio  Science, 
Vol.  12,  No.  1,  pp.  75-87,  January-February  1977. 

2 Pappert,  R.A.,  and  Goodhart,  C.L.,  "A  Numerical  Study  of  Tropospheric 
Ducting  at  HF,"  Radio  Science,  Vol.  14,  No.  5,  pp.  803-81  3, 
September-October  1979. 
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qic/n^,  but  it  might  affect  the  value  of  the  fields.  The  effect  of  profile 

normalization  on  the  calculated  field  is  discussed  in  Section  7. 

Comparison  of  DUCT  predictions  and  measurements  are  given  at  frequencies 
of  65  MHz,  170  MHz,  520  MHz  and  3300  MHz;  at  distances  111.2  km  and  222.4  km; 

for  receiver  heights  30.5  m  and  152.4  m;  and  for  transmitter  heights  that 

varied  continuously  from  the  ground  to  about  500m.  Computationally,  it  is 
more  convenient  to  set  the  transmitter  heights  at  30.5  m  and  152.4  m  and  have 
the  receiver  heights  vary  continuously.  It  was  verified  that  the  computed 
results  using  these  propagation  circuits  were  identical  to  those  for  which  the 
transmitter  and  receiver  locations  are  reversed,  as  should  be  the  case  from 
Equation  68. 

The  comparison  of  DUCT  predictions  with  mesurements  for  different 
parameters  are  found  in  Figures  16  to  19.  The  figures  denote  the  location 
(height)  of  the  duct  by  a  vertical  dashed  line  and  the  height  of  the 
transmitter  by  a  larqe  asterisk.  Also  denoted  under  the  heading  "normal"  is 
the  range  of  relative  fields  over  the  heights  shown  that  would  be  obtained  if 
the  duct  were  not  present.  This  field  is  the  median  field  due  to  troposcatter 
effects  and  was  taken  from  Reference  5. 

Figures  16  to  19  show  an  overall  excellent  agreement  between  measurements 
and  predictions,  particularly  since  the  prediction  model  is  based  on  idealized 
assumptions  (e.g.,  lateral  homogeneity)  that  only  approximate  reality.  The 
predictions  are  significantly  closer  to  the  measured  data  than  they  are  to  the 
corresponding  median  troposcatter  fields.  In  Section  7,  the  divergence  of  the 
power-sum  solutions  from  the  mode  sum  solutions  for  higher  frequencies  and  at 
higher  altitudes  is  discussed. 
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Figure  16.  Comparison  of  predictions  with  measurements  for  surface  duct, 
f  =  65  MHz.  (Mode  sum  and  power  sum  predictions  are  identical 
for  the  bottom  two  figures.) 
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Figure  17.  Comparison  of  predictions  with  measurements  for  surface  duct, 
f  =  170  MHz. 
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Figure  19.  Comparison  of  predictions  with  measurements  for  surface  duct, 
f  =  3300  MHz . 
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1'ho  mode-sum  predicted  v,i  1  nos  m  Kiautes  li>  to  I  •*  aie  essentially 
identical  to  the  predicted  values  ot  I'nppert  and  iV*odhart  (References  s  and 
_V>  .  Hie  methixi  ot  Pappei  t  and  iloodharf,  however,  reunites  spec  1 1  i  oa  t  .ion  ot  a 
reteience  height  P  at  which  the  retleotion  ooet 1 10 vent s  ot  tlf  nation  1  JO  are 
calculated.  lliey  indicate  that  us  mu  their  tormiilatn'ii  or  the  problem,  their 
results  oan  he  ditterent  tor  different  reference  heights.  In  Reference  s, 
they  present  two  examples  m  whioh  the  promotions  us  mo  p»  0  ate  ditterent 
from  those  us  mo  P  -  UO.s  m.  Ihese  are  shown  in  t'i  ante  .'0  itt  whioh  they  ate 
compared  with  the  results  us  mu  the  Pik'T  prooram.  In  both  oases ,  the  PIK'T 
predictions  are  m  aureemettt  with  the  p  -  IS.'.  4  m  results  whioh  IHpj'ert  and 
doodhart  state  are  preterahle  to  the  P  -  P  results  .nut  whioh  autee  with  the 
measurements  more  closely. 

Fd.KV.Yl’KP  Pl’i'TS 

Measurements  m  an  environment  oontammu  an  elevated  duct  were 
documented  hy  Ski  liman  and  Woivls  i  Reference  41.  lliey  compared  their 
measurements  with  predictions  obtained  us i nu  the  method  or  tttppert  and 
ii'oodhat  t  (Reference  M.  llie  ret  fact  in  tv  used  here  tv'  calculate  the  tields  is 
illustrated  tn  Flume  J  1 .  lire  retractivrty  is  the  same  as  that  used  by 
Sk  i  l  Imun  and  Woods  as  an  approximat ion  ot  the  profile  observed  over  the 
ptorMuation  pvith  durmu  the  period  m  which  the  measurements  were  made.  It  is 
norma 1 1 -ted  to  zero  at  the  uround. 

I'omjMi  isons  ot  predictions  with  measurements  are  utven  itt  Kiuure  tor  t 
=  44'*..’  MHz  and  in  Fiuure  l  tor  t  *  „V0I.’  MHz.  Ttte  predict  ions  obtained  tot 
the  44'>.  -MHz  case  were  the  same  as  those  obtained  by  Ski  liman  and  Wxods . 
However,  Ski  liman  and  Woixls  did  not  obtain  pred  ict  ions  tot  the 
_\'01  .  ’-MHz  ease  . 

Hie  overall  .iijreement  between  the  measuiements  and  pi  edict  ions  shown  ill 
KiUUies  and  .M  is  cons  idei  ed  excellent.  lire  divergence  ot  the  (X'woi  sum 
solution  1 1  om  the  mode  sum  lolut  n'n  m  the  lemon  above  the  duct  is  even  mot  e 
pronounced  in  these  results  tor  an  elevated  duct  than  they  were  tot  the 
stu  •  race-due  t  eases  considered  above.  Tti  i  s  is  discussed  further  m  Sect  ion  '. 
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SECTION  6 

EIGENVALUE  ANALYSIS 


GENERAL 


In  Section  3,  a  method  is  described  for  determining  the  roots,  or 
eigenvalues,  of  the  modal  Equation  91  in  a  complex  space,  referred  to  as  q-jg- 
space.  In  this  section,  it  will  be  seen  that  the  locus  of  the  eigenvalues  in 
q^-space  have  certain  distinct  characteristics  that  can  be  predicted  from  the 
refractivity  profile.  The  roots  in  different  portions  of  this  locus  may  be 
associated  with  different  types  of  contributions  to  the  fields  in  a  duct. 

The  correspondence  between  the  locations  of  roots  in  9-space  (see 
Section  3)  and  specific  types  of  field  contributions  (i.e.,  trapped  waves, 
chordal  waves,  multi-hop  waves)  was  discussed  by  Skillman  and  Woods  (Reference 
9).  As  will  been  seen  below,  use  of  the  q1Q-space  for  this  purpose  along  with 
the  formulation  of  the  problem  utilized  in  Section  2  make  possible  the 
prediction  of  the  location  of  the  roots  that  contribute  to  trapped  waves, 
leaky  (chordal)  waves  and  multi-hop  waves. 

Tb  determine  the  field  contribution  of  a  given  (say  nth)  eigenvalue  to 
the  field,  the  corresponding  term  in  the  field  sum  in  Equation  66  is  used. 
Thus,  the  field  contribution  of  the  nth  eigenvalue  is  given  by: 


A 

n 


=  6 

o 


X  E  e 
n  n 


(1  95) 


where  the  definitions  of  the  terms  are  given  in  Section  2.  From  the 
definition  of  the  "power  sum",  Equation  67b,  it  is  seen  that  the  contribution 
of  the  nth  eigenvalue  to  the  "power  sum"  field  is  identical  to  Equation  195. 
Thus,  separate  data  for  "mode"  and  "power"  results  are  not  needed. 
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The  parameter  En  is  a  function  of  the  receiver  height  zR  and  the 
transmitter  height  z^,  and  can  in  principle  be  cast  in  the  form  of 
Equation  68: 


E 

n 


u  ( z  )  u 
n  R  n 


(ZT> 


(1  96) 


In  order  to  determine  the  contribution  of  the  nth  mode  to  the  field  at  a 
receiver  height  zR,  a  transmitter  height  zT  must  be  specified.  It  is  seen 
from  Bguation  196  that,  regardless  of  the  value  of  zT,  the  variation  of  En 
(and  therefore  of  A^)  will  be  the  same  except  for  a  factor  un(zT).  When 
expressing  the  result  in  dB,  the  variation  of  ADB  as  a  function  of  zR  will  be 
the  same  for  different  values  of  z T,  except  for  an  additive  constant. 

COMPUTED  EIGENVALUE  LOCATIONS  -  ELEVATED  DUCTS 


To  determine  the  locations  of  the  eigenvalues  in  the  q1Q-plane,  the 
refractivity  profile  shown  in  Figure  24  will  be  assumed.  The  eigenvalues  are 
determined  using  the  methods  described  in  Section  3,  for  four  different 
propagation  frequencies:  149  MHz  (Figure  25),  449  MHz  (Figure  26),  2.2017  GHz 
(Figure  27)  and  10  GHz  (Figure  28).  (It  is  emphasized  here  that  the 
eigenvalues  are  independent  of  the  location  of  the  transmitter  and  the 
receiver.)  In  each  of  these  figures,  the  location  of  an  eigenvalue  is 
indicated  by  a  dot  on  a  graph  of  the  q^-plane  in  a  portion  of  the  figure. 

When  the  dots  are  sufficiently  close  together  on  the  scale  used,  they  will 
give  the  appearance  of  a  continuous  line  or  area.  The  eigenvalue  plot  in 
Figure  27  was  shown  in  greater  detail  in  Figure  11. 

In  each  of  the  eigenvalue  plots,  the  following  regions  may  be  identified: 

1.  A  region  in  which  the  eigenvalues  lie  on  or  near  the  real  axis 
( region  A  in  Figure  11). 
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2.  A  region  in  which  the  form  of  the  locus  of  the  eigenvalues  is  "wavy," 
but  in  which  the  overall  trend  is  diagonal  from  lower  left  to  upper  right 
(region  D  in  Figure  11).  The  lower  left  origin  of  this  region  corresponds 
approximately  to  Re(q1Q)  =  0. 

The  locus  of  the  eigenvalues  connecting  the  two  regions  above  may  be 
identified  as  consisting  of  two  additional  regions:  one  which  has  a  wavy  form 
and  extends  from  Refq-jg)  =  0  in  the  negative  real  direction  (denoted  as  region  C 
in  Figure  11);  and  finally  a  region  (denoted  as  region  B  in  Figure  11) 
connecting  the  eigenvalues  on  the  real  axis  with  region  C.  The  locus  of  the 
eigenvalues  in  region  B  is  approximately  linear. 


FIELD  CONTRIBUTIONS  FROM  EIGENVALUES  IN  DIFFERENT  REGIONS  OF  THE  q1Q-PLANE 


For  the  refractivity  profile  and  frequencies  considered  above,  it  was 
shown  that  the  locus  of  the  eigenvalues  in  the  q1Q-plane  is  always  composed  of 
well-defined  sections,  or  regions.  It  will  now  be  shown  that  all  eigenvalues 
in  the  same  region  will  contribute  to  the  propagated  fields  in  a  similar 
manner.  This  is  accomplished  in  Figures  25  to  28  by  indicating  the  field 
contribution  of  various  eigenvalues.  In  each  of  these  figures,  the  graph  of 
the  field  contribution  versus  height  corresponding  to  a  given  eigenvalue  is 
joined  to  that  eigenvalue  (in  the  q^-plane)  by  an  arrow. 


f=1  49  MHz 


The  149-MHz  case  is  considered  first.  The  four  different  regions  of  the 
locus  of  the  eigenvalues  (Figure  25)  are  identifiable  in  this  case,  but  they 
are  not  as  distinct  as  for  the  higher  frequency  cases.  The  field 
contributions  for  modes  in  each  of  these  regions  are  illustrated  as  a  function 
of  height  in  the  figure. 

In  Figure  25,  the  eigenvalue  that  lies  near  the  real  q^-axis  and 
corresponds  to  plot  A  contributes  principally  within  the  duct  with  a  small 
amount  of  leakage  above  the  duct  and  even  less  below  it.  The  single  "bump"  in 
the  curve  in  plot  A  indicates  that  the  corresponding  eigenvalue  represents  the 
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"fundamental"  waveguide  mode  of  propagation. 

The  eigenvalue  corresponding  to  plot  B  in  Figure  25,  which  lies  in  what 
has  been  identified  as  "region  B,"  produces  a  field  which  is  as  strong  above 
as  it  is  within  the  duct.  It  falls  off  quickly,  however,  as  the  receiver 
height  decreases.  The  number  of  "bumps"  in  the  curve  indicates  that  the 
corresponding  eigenvalue  represents  a  waveguide  mode  of  fourth  order. 

Plot  C  in  Figure  25  shows  the  field  contribution  due  to  an  eigenvalue  in 
region  C.  Here,  the  field  is  largest  above  the  duct  and  falls  off  slowly  as 
the  receiver  height  is  decreased  until  a  point  is  reached  at  which  the  field 
decreases  much  more  rapidly. 

The  effect  of  an  eigenvalue  in  region  D  is  shown  in  plot  D,  which  is 
similar  to  plot  C  except  for  the  fact  that  in  plot  D,  the  decrease  in  field 
strength  with  a  decrease  in  altitude  is  approximately  the  same  at  all  heights, 
whereas  in  plot  C  there  is  a  sudden  drop-off  prior  to  reaching  the  ground. 

f =449  MHz 


For  the  449  MHz  case,  the  field  contributions  are  illustrated  in 
Figure  26  (plots  A^ '  ' ,  B,  C,  D)  for  six  eigenvalues.  It  is  seen 

that  the  eigenvalues  contributing  to  plots  A^ 1  ^ ^  ,  and  A^'  represent  the 
fundamental,  second  order  and  third  order  waveguide  modes,  respectively.  In 
each  of  these  cases,  the  fields  due  to  these  eigenvalues  are  negligible 
everywhere  except  within  the  duct. 

The  field  in  plot  B  is  shown  to  be  largest  above  the  duct,  but  falls  off 
very  rapidly  at  a  height  just  below  the  duct  bottom.  The  field  shown  in  plot 
C  is  similar  to  that  in  plot  B,  but  the  “leakage”  reaches  further  below  the 
duct  before  it  falls  off  rapidly.  The  rapii  fall-off  disappears  entirely  in 
the  field  of  plot  D  in  which  the  ground  presence  appears  to  play  a  role. 
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f =2201 . 7  MHz 


The  fields  shown  in  plots  A^ 1  ^ ,  A^  2  ^ ,  A^  3  ^  and  A^4^  in  Figure  27  are 
negligible  everywhere  except  within  the  duct.  Plots  A^1^,  A^2^  and  A  ^ 3  ^ 
represent  the  lowest  order  waveguide  modes. 

The  fields  shown  in  plots  B^,  B  ^ 2  ^  and  3  ^  are  "leaky"  above  the  duct, 
but  little  or  no  leakage  occurs  below  the  duct.  Comparing  plots  1  ^ 2  ^ , 
C^3^  and  C(4)  with  the  locations  of  their  corresponding  eigenvalues,  it  is 
seen  that,  as  Re(qin)  increases,  the  height  decreases  at  which  the  respective 
fields  fall  off  rapidly.  Once  the  eigenvalue  region  corresponding  to  plots 
dO),  d(2),  d(3)  and  d(4)  ig  reached,  there  is  no  rapid  fall-off  of  the  fields 
at  any  height.  This  is  interpreted  as  indicating  that  the  ground  contributes 
to  these  fields  through  the  mechanism  of  reflection. 

f=10  GHz 


The  field  behavior  indicated  in  plots  A^ 1  ^ ,  A^  2  ^ ,  B,  1  ^ 2 ^  and  D  of 
Figure  28  are  seen  to  be  entirely  analogous  to  those  already  observed  in  the 
field  plots  (Figure  27)  for  corresponding  eigenvalue  regions  for  the  2201 -MHz 
case . 


Based  on  observations  for  the  four  frequencies  considered,  it  is 
concluded  that,  when  an  elevated  duct  is  present,  the  portion  of  the 
eigenvalue  locus  on  or  near  the  real  axis  contains  the  eigenvalues  that 
describe  "trapped  waves"  existing  almost  exclusively  within  the  duct.  _.pSince 
the  modal  attenuation  with  distance  is  given  from  Equation  195  as  |e  n  |  , 


this  attenuation  depends  on  Imfp^) 


From  Equation  160  it  is  seen  that 
(n) 


I m(p  )  would  be  very  small  for  values  of  q,n'  '  on  or  near  the  real  axis. 

n  i  o 

trapped  waves  within  the  duct,  therefore,  propagate  with  little  or  no 
attenuation  per  unit  length. 


The 


As  the  location  of  the  eigenvalues  in  the  q^g-plane  moves  in  the  positive 
real  direction,  the  eigenvalues  begin  to  represent  leaky  waves,  first  to  the 
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region  above  the  duct  and  then  to  the  region  below  the  duct.  The  leaky  waves 
below  the  duct  correspond  to  the  chordal  waves  referred  to  by  Ski liman  and 
Woods  (Reference  9).  The  leakage  reaches  lower  altitudes  as  Reiq^)  increases 
until,  in  the  neighborhood  of  Re(q1Q)  =  0,  the  leakage  reaches  the  ground. 
Beyond  this  point,  ground  reflection  becomes  important,  and  the  modes  can  be 
associated  with  "multi-hop"  waves. 

The  mathematical  basis  for  Re(q1Q)  =  0  to  represent  the  boundary  in 
q.|Q-space  between  eigenvalues  which  produce  ground-influenced  fields  (multi¬ 
hop  modes)  and  eigenvalues  which  produce  ground-independent  fields  (chordal 
modes)  is  presented  below.  Also  presented  is  a  criterion  for  determining  the 
limits  of  Re(q1Q)  in  the  q^-space  enclosing  the  region  of  the  eigenvalue 
locus  representing  trapped  modes  (region  A). 

BOUNDS  ON  TRAPPED  AND  MULTI -HOP  EIGENVALUES 


The  eigenvalues  of  interest  will  generally  lie  in  a  region  of  q1Q-spaee 
of  the  form  shown  in  Figure  29,  where  the  real  dimension  is  much  greater  than 
the  imaginary  dimension.  This  was  seen  to  be  the  situation  in  the  cases 
considered  above.  Also  shown  in  Figure  29  is  the  line  representing 
arg(q1Q)  =  2ir/3  •  As  discussed  in  Section  3,  the  eigenvalues  q10^nl  are  4:116 
zeroes  of  the  function: 


GiltV  ' 


S<q10' 


l11*  421 ' 


q22'*’ 


*qL,L-1 * 


(1  97) 


where  G1  is  the  value  of  the  determinant  in  Equation  136.  The  q —  are  related 
to  q i q  through  Equations  95,  96  and  155: 


qij 


qioti 


+  s .  . 
13 


1 16 


(198) 
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(199) 


S.  .  = 
30 


2  M(2 . )  x  10  ,  S.  real 

3 


(  200) 


and  the  refractivity  profile  is  assumed  normalized  to  zero  at  the  ground. 

Each  parameter  enters  Equation  197  through  K^q^j),  m  =  1,2,  where 

K  is  a  linear  sum  of  modified  Hankel  functions.  In  APPENDIX  A  it  was  shown 
m 

that,  if  Im(q^ j )  >  0,  the  asymptotic  behavior  of  each  Km(q^j)  for 

arg(q_")  <  2ir/3  is  different  from  its  asymptotic  behavior  for 

arg(q.  )  >  27T/3.  Since  Im(q  )  >  0,  it  follows  from  Equations  198  to  200 
ij  10  — 

that  Im(q^.)  _>  0,  so  that  the  prior  statement  is  valid  for  the  case  under 
consideration.  Therefore,  it  is  reasonable  to  assume  that,  for  a  particular 
q^  ,  the  behavior  of  the  function  G^  in  Equation  197  for  arg(q_)  <  2tt/3 
would  be  different  from  the  behavior  of  G^  for  arg(q_)  >  2rr/3.  If  the 

functional  behavior  of  G^  is  different  when  crosses  the  line 

arg(q_  )  =  2tt / 3,  then  it  is  reasonable  to  assume  that  the  locus  of  the  roots 
of  Gj  would  be  different  as  well. 


Consider,  for  example,  the  refractivity  profile  shown  in  Figure  24,  with 
a  propagation  frequency  f=  2201.7  MHz.  Using  Equations  198  to  200,  Figure  30 
illustrates  the  locations  of  the  q^  for  two  values  of  q^,  such  that  only  the 
q  i  q  (and  not  the  other  q^)  are  on  different  sides  of  the  line  arg(q)  =  2n/3. 
The  behavior  of  G^  and  the  locus  of  its  roots  would  be  expected  to  be 
different  for  the  two  cases  shown  in  Figure  30.  Similarly,  the  locus  of  the 
roots  of  Gj  are  expected  to  be  different  for  the  two  cases  shown  in  Figure  31. 


Now  if,  in  the  vicinity  of  arg(q  )  =  2n/3  ,  the  roots  of  G^fq^)  occur 
for  values  of  q1Q  for  which  Sn(q10)  is  small,  then  Re(q1Q)  would  also  be  small 
for  these  roots.  Since  a  different  behavior  of  the  locus  of  the  roots  of 
G(q1Q)  would  be  expected  for  values  of  q^Q  on  either  side  of  the  line 
arg(qio)  =  2ti/ 3,  this  difference  of  behavior  would  appear  to  occur  near 
Re ( q 1 q )  =0.  This  would  explain  the  location  in  Figure  11  of  the  boundary 
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between  region  C  and  region  D  of  the  locus  of  the  roots.  A  similar  difference 
in  behavior  across  Re(q1Q)  =  0  is  observed  in  Figues  25,  26  and  28  as  well. 
Since  q1(-j  is  the  value  of  q  (Equation  26)  at  the  ground,  it  is  not  at  all 
surprising  that  the  field  contributions  of  the  roots  for  which  Refq-jg)  >  0 
would  demonstrate  a  ground  effect  (e.g.,  plot  D  of  Figure  28),  whereas  field 
contributions  of  the  roots  for  ReCq^)  <  0  would  not  demonstrate  a  ground 

i  ~>  \ 

effect  (e.g.,  plot  C  of  Figure  28). 


The  same  reasoning  would  be  expected  to  apply  to  the  fields  within  the 
duct.  As  demonstrated  in  the  previous  subsection,  these  fields  would  receive 
their  major  contribution  from  the  "trapped"  modes,  the  eigenvalues  q^^0^ 
of  which  lie  on  or  near  the  real  q^  axis.  The  duct  boundaries  in  q-space  are 
characterized  by  the  values  of  q  corresponding  to  the  points  A  and  B  in  the 
modified  refractivity  profile  of  Figure  24.  These  values  of  q  are: 


(1)  Point  As  either  q^ 1  or  q9l 

(2)  Point  B:  either  q20  or  q32 

Since  the  eigenvalues  contributing  to  the  trapped  waves  in  the  duct 
satisfy: 


lm(q  )  a  0,  trapped  modes 

then  from  Equations  198  to  200, 


(201  ) 


Im( q .  .  )  »  0, 
13 


trapped  modes 


(202) 


It  follows  that,  for  trapped  modes,  the  roots  q^^n^  are  such  that,  if  any  q^ 
crosses  the  line  arg(q)  =  2tt/3,  it  does  so  in  the  vicinity  of: 


Re ( q .  . 

13 


0, 


trapped  modes 


(203) 


It  follows  from  Equations  202  and  203  that,  for  trapped  modes,  the 

roots  are  such  that  if  any  q^  crosses  the  line  arg(q)  =  2ir/3,  it  does 

so  in  the  vicinity  of: 
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=  0 


(204) 


From  Equations  198  to  200,  the  value  of  q^Q  at  which  Equation  204  holds  is: 


2/3 


0 


tan  a 


2  M(z . )  x  10 
D 


-6 


(205) 


But  since  the  right  side  of  Equation  205  is  dependent  only  on  the  index  j  and 
not  on  the  index  i,  the  value  of  is  the  same  for  the  two  q^ ' s  that 

characterize  point  A  in  Figure  24  (i.e.,  q1 1  and  q21  ) ,  an^  is  t^ie  same  f°r  t*le 
two  q^'s  that  characterize  point  B  (i.e.,  q22  and  q32). 

The  locations  of  the  roots  q^1^  which  represent  the  fields  within  the 
duct  would  thus  be  expected  to  be  affected  by  the  value  of  q1Q  for  which  q1 1 
( or  q  21 1  =  0i  and  by  the  value  of  q1Q  for  which  q22  (or  q32)  =  0.  To 
investigate  the  relationship  between  these  values  of  q1Q  and  the  bounds  of  the 
region  in  the  q^-plane  containing  the  eigenvalues  representing  trapped  modes, 
consider  TABLE-1 .  The  information  in  the  table  is  for  the  refractivity 
profile  and  frequencies  used  to  obtain  Figures  25  through  28.  In  the  table, 
trapped  modes  are  taken  as  those  modes  for  which  Im(q  )  <  .02. 

It  is  seen  from  TABLE-1  that,  for  the  elevated  duct  profile  and  for  the 
frequencies  considered,  the  left  boundary  4n  the  q^  plane  of  the  q^11^ 
representing  trapped  waves  in  the  duct  is  approximated  well  by  the  value  of 
q10  for  which  q^  =0.  Similarly,  the  right  boundary  of  these  trapped  modes 
is  approximated  well  by  the  value  of  q1Q  for  which  q22  =  0. 

It  is  also  interesting  to  note  from  TABLE-1  that  the  number  of  trapped 
modes  between  the  right  and  left  boundaries  of  the  trapped-wave  eigenvalues 
(density  of  the  eigenvalues  in  q^-space)  increases  only  logarithmically  with 
frequency.  Therefore,  mesh  size  in  the  numerical  method  descibed  in  Section  3 
to  locate  the  roots  need  not  be  altered  appreciably  for  different  frequencies. 
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APPLICATION  OF  EIGENVALUE  ANALYSIS  TO  TOO-DUCT  PROFILE 


In  order  to  verify  the  concepts  described  above  for  estimating  the 
locations  of  trapped  modes,  the  theoretical  two-duct  geometry  illustrated  in 
Figure  32  will  be  considered  with  a  propagation  frequency  of  449  MHz.  The 
computed  eigenvalues  are  illustrated  as  the  dots  in  the  q 1(^  plane  in  Figure 
33.  The  eigenvalues  in  this  figure  appear  to  lie  along  two  separate  loci,  one 
characterizing  each  duct.  The  trapped  waves  for  each  duct  are  easily 
identifiable  as  those  lying  near  the  real  axis.  Field  plots  in  Figure  33 
illustrate  the  field  contributions  for  each  of  the  labelled  different 
eigenvalues.  (Note  that  the  labelling  is  not  in  any  apparent  sequence,  but 
rather  represents  the  number  of  the  mode  in  the  order  in  which  it  is 
calculated  using  the  MODESRCH  method  described  in  Section  3) . 

In  Figure  33,  plots  82,  81  and  80  are  shown  to  correspond  to  the  trapped 
waves  in  the  higher  duct.  Plot  56  represents  a  leaky  wave  from  the  same 
duct.  Plots  44,  46,  and  47  are  shown  to  correspond  to  trapped  waves  in  the 

lower  duct.  It  is  interesting  to  note  that  the  eigenvalue  corresponding  to 
plot  82  (i.e.,  the  leftmost  eigenvalue)  represents  the  fundamental  waveguide 
mode  for  the  system,  and  contributes  to  the  field  in  the  upper  duct.  Plot  44, 
on  the  other  hand,  represents  a  higher  mode  of  the  system,  but  is  the 
fundamental  mode  of  propagation  in  the  lower  duct. 


Using  Equation  205,  the  boundaries  of  the  region  in  the  q1Q  plane 
containing  the  trapped  modes  of  the  lower  duct  are  approximated  by: 


1  1 


‘10 


=  -  11.7 


and 


22 


l10 


=  -  7.  3 


1  24 
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while  the  corresponding  actual  boundaries  for  the  eigenvalues  shown  in 
Figure  33  are  -10.8  and  -8.1.  The  boundaries  of  the  region  containing  the 
trapped  modes  of  the  upper  duct  are  approximated  by: 


q 


10 


q 


10 


1  8.  4 

1  4.  1 


while  the  corresponding  actual  boundaries  for  the  eigenvalues  are  -1 7. 7  and 
-1  4.  5. 


The  method  described  in  the  preceding  subsections,  for  approximating  the 
boundaries  of  the  region  in  the  q1(^  plane  that  contains  trapped  modes,  is  thus 
seen  to  have  a  more  general  validity.  Its  application  to  ground-based  ducts 
will  be  demonstrated  in  a  later  subsection. 

COMPUTED  EIGENVALUE  LOCATIONS  -  SURFACE  DUCTS 


The  locations  of  the  eigenvalues  in  the  q1Q-plane  will  now  be 
demonstrated  for  the  surface  duct  refractivity  profile  illustrated  in 
Figure  34  at  four  propagation  frequencies:  65,  170,  520,  and  3300  MHz. 

The  eigenvalue  locations  for  the  u  1Hz  case  is  illustrated  in  Figure 
35.  Unlike  the  elevated  cases,  the  locus  of  the  eigenvalues  appear  to  have 
two  branches:  one  including  the  eigenvalue  which  produces  plot  C  and  one  not 
including  this  eigenvalue.  The  field  plots  shown  in  Figure  35  illustrate  the 
field  contributions  of  each  of  the  modes.  The  modes  that  appear  to  be 
affected  by  the  presence  of  the  duct  are  those  corresponding  to  plots  A  and 
C.  Thus,  in  Figure  35,  it  would  be  reasonable  to  assume  that  plots  A  and  C 
correspond  to  eigenvalues  that  belong  to  another  branch.  Since  the  former 
locus  is  associated  with  the  ducted  fields,  it  would  be  expected  that  the 
latter  locus  is  produced  by  a  different  propagation  mechanism.  It  will  be 
demonstrated  below  for  the  170-MHz  case  that  this  mechanism  is  likely  to  be 
simple  diffraction. 
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The  eigenvalue  locations  for  the  170-MHz  case  are  illustrated  in 
Figure  36.  Again,  two  distinct  loci  may  be  identified:  one  consisting  of 
eigenvalues  corresponding  to  plots  A,  B,  C,  D,  E,  F,  G;  and  the  other 
consisting  of  eigenvalues  corresponding  to  plots  H,  I  and  J.  Hie  fields  of 
plots  A,  B,  C,  D,  E,  F  and  G  are  seen  to  be  affected  by  the  duct.  However, 
the  fields  in  plots  H,  I,  and  J,  seem  unaffected. 

It  is  interesting  to  note  from  Figures  35  and  36  that  the  loci  of  the 
modes  affected  by  the  duct  have  a  positive  slope,  while  the  loci  of  the  modes 
not  affected  by  the  duct  have  a  negative  slope.  Tb  ascertain  the  likely 
source  of  the  fields  that  are  not  affected  by  the  duct,  the  "ductless" 
refractivity  profile  shown  in  Figure  37  will  be  considered.  The  eigenvalues 
produced  by  this  profile  for  f  =  1 70  MHz  are  shown  in  Figure  38.  Hie  form  of 
the  field  contribution  for  each  of  these  eigenvalues  is  the  same,  a  typical 
one  being  illustrated  by  a  single  field  plot  in  the  figure.  Hie  form  of  this 
field  is  seen  to  resemble  that  of  those  field  plots  in  Figures  35  and  36  which 
do  not  appear  affected  by  the  presence  of  the  duct.  Hie  slope  of  the  locus  of 
eigenvalues  in  Figure  38  is  seen  to  be  negative,  resembling  the  loci  of  the 
eigenvalues  that  are  not  affected  by  the  duct  in  Figures  35  and  36.  But  when 
the  refractivity  profile  does  not  exhibit  a  duct,  the  fields  at  beyond-horizon 
distances  are  caused  by  diffraction  (since  troposcatter  has  been  ignored  in 
this  formulation) .  The  eigenvalues  in  Figures  35  and  36  that  appear  to  be 
transparent  to  the  duct  are,  therefore,  most  likely  associated  with 
diffraction  fields. 

Resuming  the  investigation  of  the  eigenvalue  locations  for  surface  duct 
environments,  Figure  39  illustrates  such  locations  for  a  propagation  frequency 
of  520  MHz  and  the  duct  profile  shown  in  Figure  34.  One  eigenvalue 
(corresponding  to  plot  D)  is  now  observed  which  does  not  lie  on  the  locus  of 
the  other  eigenvalues.  Hie  eigenvalues  corresponding  to  plots  B  and  C,  which 
do  not  lie  near  the  real  q1Q  axis,  are  seen  to  be  highly  leaky  to  the  region 
above  the  duct.  Hie  field  in  plot  D,  which  corresponds  to  an  eigenvalue  that 
does  not  lie  on  the  locus  of  the  other  modes,  has  the  appearance  of  the 
diffracted  field  shown  in  Figure  37  except  for  an  oscillation  in  the  curve 
near  the  ground. 

131 


Figure  39.  Eigenvalues  and  some  of  their  field  contributions  for  refractivity  profile 
of  Figure  34,  f  =  520  MHz. 


ESD-TR-81 -1 02 


Section  6 


The  eigenvalues  calculated  for  a  frequency  of  3300  MHz  are  illustrated  in 
Figure  40.  The  field  contributions  of  the  first  three  trapped  modes  are  shown 
in  plots  A^,  A^  and  A^3'.  The  contribution  of  a  mode  located  off  the  real 
axis  is  illustrated  in  plot  B.  It  is  seen  that  this  mode  still  contributes 
significantly  to  the  field  within  the  duct,  and  leaks  slightly  above  the  duct 
and  for  only  a  short  distance. 

Just  as  in  the  case  of  elevated  ducts  (TABLE  1),  the  prediction  method 
(Equation  205)  for  approximating  the  bounds  of  Refq^)  for  the  trapped  modes 
in  the  q^g-plane  should  apply  for  surface  ducts  as  well.  The  predictions 
using  this  approximation  are  compared  in  TABLE  2  with  the  actual  results 
calculated  by  the  method  of  Sections  2  and  3  (and  illustrated  in  Figures  35, 

36,  39  and  40). 

Notice  that  TABLE  2  shows  no  trapped  modes  for  the  65-MHz  case,  whereas 
Figure  35  indicates  the  presence  of  a  strong  ducted  field.  There  is  no 
contradiction,  however,  since  the  trapped  modes  as  defined  for  the  purposes  of 
TABLE  2  are  non-leaky  and  are  essentially  unattenuated  along  the  waveguide. 

The  field  contribution  evident  in  Figure  35  is  leaky  and  is  attenuated  along 
the  waveguide.  It  is,  however,  appreciable  at  the  distance  for  which  the 
calculations  were  made  for  Figure  35. 


Finally,  note  that  no  extreme  change  occurs  in  the  locus  of  the 
eigenvalues  in  Figures  36,  39,  and  40  in  the  vicinity  of  Re  (q-|g)  =  0,  as  was 
the  case  for  elevated  ducts.  For  elevated  ducts,  it  was  shown  that  such  a 
change  occurs  between  modes  that  produce  fields  affected  by  the  ground,  and 
modes  that  produce  fields  unaffected  by  the  ground.  However,  for  surface 
ducts,  the  ground  plays  a  role  in  the  duct  itself  and,  therefore,  affects  the 
trapped  waves.  This  is  demonstrated  by  the  fact  that  the  trapped  modes  lie  on 


either  side  of 


0 
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INFLUENCE  OF  DUCT  HEIGHT 


To  illustrate  the  effects  of  the  duct  height  on  the  locations  of  the 
roots  of  the  modal  equation  in  q^g-space,  consider  the  refractivity  profile 
shown  in  Figure  41.  This  profile  (in  M-units)  is  identical  to  that  of  Figure 
24  except  the  center  of  the  duct  in  Figure  41  is  at  about  half  the  height  of 
the  one  shown  in  Figure  24.  The  eigenvalues  in  q^g-space  for  a  frequency  of 
449  MHz  using  the  profile  of  Figure  41  are  illustrated  in  Figure  42.  The 
corresponding  eigenvalues  using  the  profile  of  Figure  24  are  illustrated  in 
Figure  26. 

Comparing  Figures  26  and  42,  notice  that  the  number  of  trapped  modes 
(those  near  the  real  axis)  is  the  same.  These  trapped  modes  appear  at  lower 
negative  values  of  Re(q1g)  in  Figure  42  than  in  Figure  26  since,  from  Equation 
205,  the  values  of  q1Q  at  which  q^  and  q^2  vanish  have  smaller  negative 
values  in  the  profile  in  which  the  duct  height  is  lower.  Although  the  number 
of  trapped  modes  is  the  same,  the  number  of  all  other  modes  (or  the  density  of 
the  other  modes  in  q1Q-space)  is  significantly  smaller  for  the  lower  duct 
height  than  for  the  higher  duct  height.  In  addition,  there  are  fewer  and  less 
intense  fluctuations  in  the  eigenvalue  locus  of  Figure  42  than  in  the  one  in 
Figure  26.  It  is  reasonable  to  assume  that,  when  the  duct  height  decreases  to 
a  point  where  the  elevated  duct  becomes  a  surface  duct,  these  fluctuations 
would  disappear  entirely,  as  was  seen  to  be  the  case  in,  say,  Figure  40. 
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INFLUENCE  OF  DUCT  THICKNESS  AND  INTENSITY 


The  propagation  frequency  of  449  MHz  will  be  utilized  to  investigate  the 
effect  on  the  eigenvalue  locus  of  a  variation  of  duct  thickness  and  intensity. 
Figures  43  through  48  provide  illustrations  of: 

( 1 )  eigenvalue  locations  and 

(2)  modal  contributions 

for  ducts  of  decreasing  thickness  and  intensity.  In  each  case,  the  modified 
refractivity  gradients  in  each  region  are  the  same  as  those  in  the  profile  of 
Figure  24.  The  boundary  between  region  2  and  region  3,  however,  is  lowered 
until,  in  Figure  48,  region  2  disappears  entirely. 

As  would  be  expected,  Figures  43  and  48  show  that  the  number  of  trapped 
modes  decreases  as  the  duct  becomes  less  intense  and  narrower.  The 
fluctuation  in  the  eigenvalue  locus  also  becomes  subdued  and  disappears 
entirely  when  the  duct  disappears.  The  highest  order  modes  become 
increasingly  leaky  as  the  duct  strenqth  diminishes.  The  fluctuation  in  the 
modal  contributions  to  the  field  also  becomes  small,  so  that,  in  plot  ( B)  of 
Figure  48,  the  inflection  in  the  curve  (see  arrow  in  figure) ,  which  identifies 
the  contribution  as  originating  from  the  second  order  mode,  is  hardly 
distinguishable.  In  Figure  48,  the  lowest  order  mode  from  which  plot  (A)  was 
constructed  was  the  mode  for  which  Re(q1Q)  is  minimum  —  this  even  though  it 
does  not  have  the  lowest  attenuation  per  kilometer  (since  it  does  not  have  the 
lowest  value  of  Im(q1Q)).  The  contribution  of  the  mode  with  lowest 
attentuation  per  kilometer  is  shown  in  plot  (C)  of  the  figure. 


Figure  43.  Eigenvalues  and  two  of  their  field  contributions  for  a  refractivity  profile  with  the  same 
refractivity  gradients  and  duct  height  as  the  profile  of  Figure  24,  but  with  duct 
thickness  =  0.4  km. 


XMTR  HEIGHT 


45.  Eigenvalues  and  two  of  their  field  contributions  for  a  refractivity  profile  with  the 
same  refractivity  gradients  and  duct  height  as  the  profile  of  Figure  24,  but  with 
duct  thickness  =  0.2  km. 
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Figure  47.  Eigenvalues  and  two  of  their  field  contributions  for  a 

refractivity  profile  with  the  same  refractivity  gradients 
and  duct  height  as  the  profile  of  Figure  24,  but  with  duct 
thickness  =  0.05  km. 


Figure  48.  Eigenvalues  and  some  of  their  field  contributions  for  a  refractivity  profile  with  the  same 
refractivity  gradients  and  duct  height  as  the  profile  of  Figure  24,  but  duct  thickness  = 
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SECTION  7 

ANALYSIS  OF  DUCTED  FIELDS 


FIELD  SIMULATION  BY  PARTIAL  MODE  SUM 


In  the  previous  section,  it  was  shown  that  specific  types  of  propagation 
may  be  associated  with  different  eigenvalues.  Thus  some  eigenvalues  will 
contribute  only  to  trapped  waves  within  the  duct,  some  will  contribute  to 
leaky  (or  chordal)  waves  and  some  will  be  associated  with  multi-hop 
contributions  caused  by  reflections  from  the  ground.  It  was  also  shown  that 
the  locations  of  the  eigenvalues  in  q^-space  that  provide  each  type  of 
contribution  can  generally  be  identified. 

From  Equations  67a  and  68,  the  relative  field  may  be  written  as: 

-jp  r 
n 

A  =  B  |  E  X  u  (z)u  (z)e  |  (206) 

o  n  n  n  R  n  T 

where  zR  and  zT  are  the  receiver  and  transmitter  heights,  respectively,  and 

the  sum  is  over  the  eigenvalues.  It  will  be  noticed  from  Equation  206  that 

for  each  mode,  there  is  reciprocity  between  the  transmitter  and  receiver  — 
that  is,  the  result  remains  unchanged  if  the  transmitter  and  receiver  heights 
are  interchanged.  Also,  u^Zp)  characterizes  the  influence  of  the  receiver 
height  to  the  modal  field  contribution,  and  un(z^,)  characterizes  the  influence 
of  the  transmitter  height  to  the  modal  contribution.  This  implies  that,  if 
both  zR  and  zT  are  heights  at  which  the  same  type  of  wave  is  dominant  (e.g. 
trapped  wave,  chordal  wave,  multip-hop  wave),  then  the  dominant  terms  in  the 
sum  of  Equation  206  will  be  those  of  the  eigenvalues  that  contribute  most  to 

that  type  of  wave.  Thus,  if  both  zT  and  zR  are  heights  within  the  duct,  the 

terms  in  Equation  206  representing  the  eigenvalues  in  region  A  (see  Section  6) 
would  contribute  most  significantly  to  the  overall  sum.  Similarly,  if  zT  and 
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z^  are  both  near  the  ground,  the  dominant  contributions  in  Equation  206  will 
be  derived  from  the  eigenvalues  in  region  D. 

If  is  a  height  at  which  a  particular  type  of  wave  is  dominant,  and  zR 

is  a  height  at  which  a  different  type  of  wave  is  dominant,  then  the 
significant  terms  in  the  sum  in  Equation  206  would  be  those  for  the 
eigenvalues  characterizing  each  of  these  types  of  waves. 

The  above  considerations  imply  that,  for  certain  specific  propagation 
circuits  of  interest,  it  is  possible  to  describe  the  relative  fields  by  a 
relatively  small  number  of  eigenvalues.  Since  it  is  often  possible  to 
localize  the  region  of  the  q^-plane  in  which  these  eigenvalues  lie  (using  the 
methods  described  in  the  previous  sections),  a  significant  saving  in 
computation  time  may  be  realized. 

Illustrating  the  point  made  above,  consider  the  propagation  in  the  duct 
in  Figure  24.  At  2201.7  MHz,  the  resulting  eigenvalues  were  illustrated  in 
Figure  27.  A  portion  of  these  eigenvalues  is  shown  in  Figure  49.  The 
location  of  the  eigenvalues  is  independent  of  the  transmitter  and  receiver 
heights.  Assume  the  receiver  and  transmitter  are  both  located  within  the 
duct.  Then  the  dominant  contribution  to  the  field  will  derive  from  the 
eigenvalues  close  to  the  real  axis  in  Figure  49.  Figure  50  compares  the  field 
using  all  modes  for  which  ReCq^g)  <  -77  with  the  total  field  using  all  the 
modes  shown  in  Figure  27.  It  is  seen  that  within  the  region  of  interest 
(i.e.,  within  the  duct),  there  is  good  agreement  between  the  two  cases. 
However,  outside  the  duct,  a  large  discrepancy  between  the  two  results  is 
apparent,  since  eigenvalues  that  contribute  to  the  regions  outside  the  duct 
were  not  included  in  the  partial  mode  sum.  Figure  51  compares  the  field  using 
all  modes  for  which  Re(q1(j)  <  -75  (i.e.,  more  modes  than  were  used  in  Figure 
50)  with  the  total  field  using  all  modes.  The  agreement  within  the  duct  is 
even  better  than  in  the  former  case,  and  the  field  above  the  duct  using  the 
partial  sum  is  closer  to  the  total  field  than  is  the  partial  sum  result  in 
Figure  50. 
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Figure  49.  A  portion  of  the  eigenvalues  shown  in  Figure  27. 
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DIFFERENCES  BETWEEN  MODE  AND  POWER  SUMS  -  EFFECT  OF  PROFILE  NORMALIZATION 

In  Section  5  it  was  seen  that,  in  some  cases,  the  "power  sum”  result 
(Equation  67d  in  Section  2)  differed  greatly  from  the  "mode  sum"  result 
(Equation  67c).  This  difference  was  most  obvious  in  Figures  22  and  23  in  the 
region  above  the  duct,  and  was  apparently  caused  by  extreme  destructive 
interference  between  individual  modal  contributions  that  are  significant  at 
those  heights.  From  the  measured  data  shown  in  these  figures,  it  appears  that 
the  power  sum  result  is  more  realistic  than  the  mode  sum  result.  Since  the 
mode  sum  is  theoretically  more  accurate  than  the  "power  sum,  the  question 
arises  as  to  the  reason  for  the  discrepancy  between  the  measurements  and  the 
mode  sum  results. 

Skillman  and  Woods  (Reference  9)  suggested  that  the  source  of  the  problem 
might  be  a  perturbation  caused  by  horizontal  refractivity  inhomogeneities,  or 
by  inadequacy  of  the  mathematical  model.  The  fact  that  the  problem  seems  to 
occur  only  at  high  altitudes  might  point  more  to  the  latter  explanation. 
Pekeris  (Reference  4)  showed  that  the  "earth-flattening"  approximation  used  in 
this  mathematical  model  (in  which  the  curvature  of  the  earth  was  taken  into 
account  by  modifying  the  refractivity  profile)  would  become  less  accurate  as 
the  altitude  increases.  This  inaccuracy  would  affect  both  modal  amplitudes 
and  modal  phases.  Since  phases  would  generally  have  a  more  dramatic  influence 
on  the  calculated  field  than  small  differences  in  amplitudes,  the  solution 
that  strictly  accounts  for  these  phases  would  be  expected  to  go  awry  prior  to 
the  solution  that  assumes  phase  incoherency.  This  would  explain  the  reason 
for  the  power  sum  solution  providing  predictions  in  Figures  21  and  22  that  are 
superior  to  the  mode  sum  predictions. 

Tb  illustrate  the  effect  on  the  relative  field  strength  of  small 
differences  in  modal  phases,  it  is  interesting  to  compare  the  relative  field 
strength  results  for  the  case  of  a  "normalized"  refractivity  profile  (see 
Section  3),  with  the  case  of  an  unnormalized  refractivity  profile.  As 
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discussed  in  Section  3,  the  principal  computational  difference  between  these 
two  cases  appears  in  the  value  of  pn  in  the  exponent  in  Equation  206.  From 
Equation  94,  for  the  unnormalized  profile: 


p  =  k 
n  o 


( n) 


1  -  H  tan  a  - 

!  1  i 


1  0 


2/3 


tan  a. 


1 


(207) 


while  for  the  normalized  profile, 


p  =  k 
n  o 


(n) 


M  0 


2/3 


tan  a 


(208) 


Since  | tan  a  |  <<  1,  the  difference  between  pn  as  defined  by  Equation  207  and 
Pn  as  defined  by  Equation  208  is  very  small.  Figure  52  provides  results  using 
the  normalized  refractivity  profile  of  Figure  24  and  using  a  corresponding 
unnormalized  profile  for  a  frequency  of  449  MHz.  In  the  unnormalized 
refractivity  profile,  the  modified  refractivity  at  the  ground  is  taken  as  300 
(instead  of  zero).  The  power-sum  results  for  both  the  normalized  and 
unnormalized  profiles  are  identical.  However,  these  profiles  produce 
different  results  using  the  mode  sum.  Note  also  that  the  results  differ  only 
in  the  region  above  the  duct,  but  are  the  same  within  and  below  the  duct. 
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EFFECT  OF  TRANSMITTER  HEIGHT  -  SINGLE  DUCT 


Elevated  Duct 


To  investigate  the  effect  of  the  transmitter  height  on  the  relative  field 
strength,  the  elevated  duct  profile  of  Figure  24  will  be  used  for  a  frequency 
of  2201.7  MHz.  The  mode-sum  and  power-sum  results  are  presented  in  Figure  53 
in  order  of  decreasing  transmitter  height  (indicated  by  an  asterisk  in  the 
left  portion  of  each  plot.)  In  Fiqures  53a  and  53b,  the  transmitter  is 
located  above  the  duct.  The  transmitter  is  within  the  duct  in  Figures  53c, 
53d,  53e ,  and  53f .  In  Figures  53g  and  53h,  the  transmitter  is  below  the  duct. 

When  the  transmitter  is  above  or  in  the  duct,  the  fields  in  or  above  the 
duct  are  relatively  large.  The  optimum  coupling  of  energy  into  the  duct 
occurs  when  the  transmitter  is  near  the  duct  center.  When  the  transmitter  is 
above  or  in  the  duct,  the  fields  are  greatest  within  the  duct.  When  the 
transmitter  is  below  the  duct,  field  enhancement  due  to  the  presence  of  the 
duct  is  negligible  or  absent  entirely. 

By  reciprocity,  the  results  cited  above  would  hold  if  the  receiver  and 
transmitter  were  interchanged.  It  may  therefore  be  concluded  that,  under 
suitable  duct  environments  (i.e.,  environments  in  which  the  minimum  trapping 
frequency  of  the  duct  is  less  than  the  transmission  frequency),  the  duct 
serves  to  significantly  enhance  the  fields  when  both  the  transmitter  and 
receiver  are  in  or  above  the  duct.  If  either  terminal  of  the  propagation 
circuit  is  below  the  duct,  the  resulting  fields  would  not  be  significantly 
enhanced  by  the  presence  of  the  duct. 

Surface  Duct 


The  fields  in  a  surface  duct  will  now  be  considered  for  different 
transmitter  heights.  The  duct  profile  of  Figure  34  will  be  used  for 
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a  frequency  of  3300  MHz.  Fiqure  54a  shows  the  fields  calculated  using  the 
mode-sum  series  for  three  transmitter  heights:  one  outside  the  duct  and  two 
within  the  duct.  Figure  54b  shows  the  corresponding  fields  using  the  power 
sum  solution.  Notice  that,  contrary  to  the  case  of  elevated  ducts,  the  field 
is  weaker  in  the  duct  than  above  it  when  the  transmitter  is  above  the  duct. 
This  would  be  due  to  the  fact  that,  when  the  source  and  observer  are  both 
above  the  duct  and  separated  by  a  distance  of  111.2  km,  they  are  within  line 
of  sight  of  each  other. 

In  Figures  54a  and  54b,  the  transmission  is  taken  as  being  vertically 
polarized.  The  results  for  this  case  were  found  to  be  identical  to  those 
using  horizontal  polarization. 
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Figure  54a.  Mode  sum  fields  calculated  using  the  refractivity  profile  of 
I  Figure  34  for  different  transmitter  heights. 
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Figure  54b.  power  sum  fields  calculated  using  the  refractivity  profile  of 
Figure  34  for  different  transmitter  heights. 
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EFFECT  OF  TRANSMITTER  HEIGHT  -  DOUBLE  DUCT 


The  effect  of  transmitter  height  on  relative  field  strength  will  now  be 
investigated  for  the  two-duct  refractivity  profile  illustrated  in  Figure  32 
for  a  frequency  of  449  MHz.  The  mode-sum  and  power-sum  results  are  presented 
in  Figures  55a  through  55d  in  order  of  decreasing  transmitter  height.  In 
Figure  55a,  the  transmitter  is  above  the  upper  duct;  in  Figure  55b,  the 
transmitter  is  within  the  upper  duct;  in  Figure  55c,  the  transmitter  is 
between  the  upper  and  lower  ducts,-  and  in  Figure  55d,  the  transmitter  is 
within  the  lower  duct.  The  effect  of  each  duct  (or  lack  of  effect)  is  obvious 
in  each  figure.  It  is  clear  that  the  field  is  qreatest  within  each  duct  when 
the  transmitter  is  located  within  it.  From  Figures  55c  and  55d,  it  is  seen 
that,  when  the  transmitter  is  below  the  upper  duct,  the  field  in  the  upper 
duct  is  similar  to  the  field  that  would  be  expected  if  the  upper  duct  were  not 
present.  The  lower  duct,  however,  serves  to  enhance  the  field  in  these  cases . 

The  conclusions  drawn  for  the  single  duct  environment  may  be  generalized 
to  the  two-duct  environment:  Each  duct  will  significantly  enhance  the  fields 
when  both  the  transmitter  and  receiver  are  in  or  above  it.  If  either  terminal 
of  the  propagation  circuit  is  below  one  duct  (or  both  ducts),  then  the 
resulting  field  will  not  be  enhanced  by  the  presence  of  that  (those)  ducts. 
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EFFECT  OF  DUCT  THICKNESS  AND  INTENSITY 


The  propagation  frequency  of  449  MHz  will  now  be  utilized  to  investigate 
the  effect  on  the  calculated  fields  of  a  variation  of  duct  thickness  and 
intensity.  In  each  case  considered,  the  modified  refractivity  gradient  in 
each  region  is  the  same  as  that  in  the  profile  of  Figure  24.  The  variation  in 
duct  thickness  and  intensity  is  obtained  by  lowering  the  boundary  between 
region  2  and  region  3  until  region  2  disappears  entirely.  This  is 
accomplished  for  three  different  transmitter  heights  in  Figures  56,  57  and  58, 
respectively.  In  Figure  56,  the  transmitter  is  located  above  the  duct,  in 
Figure  57  the  transmitter  is  located  within  the  duct,  and  in  Figure  58,  the 
transmitter  is  located  below  the  duct.  In  each  figure,  relative  field 
calculations  are  presented  for  different  duct  thicknesses  (and  intensities) . 

In  Figures  56  and  57  (in  which  the  transmitter  is  above  the  duct  and 
within  the  duct,  respectively)  it  is  observed  that,  for  the  first  two  or  three 
duct  thicknesses  considered,  the  maximum  field  within  the  duct  remains  about 
the  same.  The  fields  in  the  duct  fall  off  rapidly  for  ...nailer  ducts.  In 
Figure  58,  the  fields  are  very  similar  for  all  finite  duct  thicknesses,  but 
are  appreciably  smaller  in  the  limit  of  an  absence  of  the  duct. 

The  discrepancy  in  the  figures  between  the  mode-sum  and  the  power-sum 
results  for  the  smaller  duct  sizes  was  not  investigated. 
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Figure  56.  Fields  calculated  using  a  refractivity  profile  with  the  same 
^activity  gradients  and  duct  height  as  the  profile  of 
Figure  24,  for  different  duct  thicknesses  (transmitter 
above  the  duct) . 
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Figure  57.  Fields  calculated  using  a  refractivity  profile  with  the  same 
refractivity  gradients  and  duct  height  as  the  profile  of 
Figure  24,  for  different  duct  thicknesses  (transmitter 
with  the  duct) . 
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Figure  58.  Fields  calculated  using  a  refractivity  profile  with  the  same 
refractivity  gradients  and  duct  height  as  the  profile  of 
Figure  24,  for  different  duct  thicknesses  (transmitter 
below  the  duct) . 
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EFFECT  OF  FREQUENCY 


The  effect  of  frequency  variation  on  the  calculated  field  will  now  be 
investigated  for  the  elevated  duct  profile  of  Figure  24.  Measurements 
documented  by  Skillman  and  Woods  (Reference  9)  for  a  duct  environment  similar 
to  this  indicated  that  the  relative  field  variation  with  height  was  similar 
for  frequencies  449.2  and  2201.7  MHz.  Their  results  were  illustrated  in 
Section  5.  They  succeeded  in  predicting  the  results  for  frequencies  of  149.3 
and  449.2  MHz  using  a  computer  program  developed  by  Pappert  and  Goodhart 
(Reference  5), 

In  Figure  59,  calculated  results  are  presented  for  frequencies  149,  449, 
2201.7  and  10000  MHz.  The  results  for  the  449  and  2201.7  MHz  cases  were 
compared  in  Section  5  with  measurements  documented  by  Skillman  and  Woods 
(Reference  9).  It  is  seen  here  that  not  only  are  the  449  and  2201.7  MHz 
results  similar  to  each  other,  but  they  are  both  similar  to  the  results  for 
10000  MHz.  This  bears  out  the  suggestion  of  Skillman  and  Woods  that,  in  duct 
environments,  it  might  be  possible  to  utilize  calculations  at  lower 
frequencies  to  determine  relative  fields  at  higher  frequencies  . 
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SECTION  8 

CONCLUSIONS  AND  RECOMMENDATIONS 

CONCLUSIONS 

1.  A  reliable  computer  model,  based  on  a  waveguide  mode 
formulation,  has  been  developed  for  calculating  beyond-line-of-sight 
electromagnetic  fields  in  a  horizontally  stratified  tropospheric  duct 
environment. 

2.  The  model  is  capable  of  evaluating  fields  for  higher  elevated 
ducts  and  higher  frequencies  (i.e.,  through  SHF)  than  were  previously  feasible 
in  other  available  models,  which  had  numerical  difficulties  that  precluded 
computation  of  all  significant  eigenmodes  of  the  system. 

3.  The  types  of  numerical  difficulties  encountered  in  other  models 
were  eliminated  by  the  use  of  a  unique  mathematical  formulation  that 

a.  Assures  linear  independence  (Equations  (04  to  106),  even  in 
a  "numerical  sense",  of  the  homogeneous  form  of  the  governing  differential 
equation;  and 

b.  Provides  flexibility  for  judiciously  choosing  the 
particular  solution  (Equation  189)  to  the  inhomogeneous  form  of  the  governing 
differential  equation. 

4.  The  model  methodically  and  efficiently  determines  all 
significant  eigenvalues  (Section  3)  and  corresponding  fields  (section  4)  for 
elevated  and  surface  ducts  at  all  frequencies  through  SHF. 

5.  Criteria  have  been  developed  (Equation  205)  for  associating 
specific  types  of  field  contributions  with  eigenvalues  in  a  specific  portion 
of  the  eigenvalue  locus.  These  criteria  have  the  potential  for  increasing  the 
computational  efficiency  of  the  model  in  certain  circumstances. 
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6.  The  model  has  been  verified  by  comparing  its  predictions  with 
measurements  in  both  elevated  and  surface  duct  environments  (Section  5).  It 
is  the  only  known  computer  model  that  provides  predictions  which  agree  within 
a  few  dB  with  field  strength  measurements  performed  in  an  elevated  duct 
environment  at  frequencies  as  high  as  2201  MHz  (Figure  23). 

7.  The  model  is  deemed  to  be  a  valid  tool  for  predicting  and 
s-udying  electromagnetic  fields  in  a  tropospheric  duct  environment. 


RECOMMENDATIONS 

1.  Measured  data  on  the  probability  of  occurrence  of  ducts  of 
various  characteristics  in  various  geographical  locations  should  be  analyzed 
and  used  in  the  deterministic  model  in  a  manner  that  would  provide  statistical 
loss  values.  The  results  should  specify  the  probability  of  loss  as  a  function 
of  geographical  locations  and  seasons.  These  statistical  loss  values  should 
be  used  to  revise  the  long-term  power  fading  subroutine  used  in  ECAC 
propagation  models. 

2.  The  model  should  be  extended  to  predict  loss  values  at  distances 
within  line-of-sight . 

3.  The  model  should  be  extended  to  account  for  the  effect  of 
obstructions  in  front  of  the  antenna. 

4.  The  model  should  be  extended  to  allow  a  nonhomogeneous 
horizontal  ref ractivity . 
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APPENDIX  A 

ASYMPTOTIC  BEHAVIOR  OF  h, (q)  AND  h2(q)  AND  F(q)  IN  THE  UPPER 

HALF  OF  THE  q -PLANE 

Values  of  the  real  and  imaginary  parts  of  h-|  in  the  complex  q-plane  are 
shown  in  Figures  A-1  and  A-2,  respectively.  The  corresponding  values  for  h2 
are  given  in  Figures  A-3  and  A-4.  When  |qj  is  small,  h1  and  h2  are  both  of 
order  unity,  as  must  he  F(q),  defined  from  Equation  98  as: 


4 

3  ^ 

F(q)  =  h2  (q)  -  e  (q)  (A-1) 


| h 1 | ,  | h 2 1  and  |f|  only  become  exponentially  large  or  small  when  |qj  is 

large.  These  functions  may  therefore  be  studied  using  their  asymptotic 
expansions.  These  asymptotic  expansions  will  be  defined  according  to 
Reference  14.  It  is  convenient  to  define  the  following: 


1 


4  -5itj/1  2  m  -3m/2 

u ( q )  =  Aq  e  [1  +  E  (-j)  C  q  ]  (A-2) 

m=1  m 


1 

-  —  5 tt j /I  2  m  -3m/2 

v(q)  =  Aq  4  e  (1+E  (j)  Cq  ) 

m=1  m 


(A-3) 
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Appendix 


Im  (q) 


Figure  A-2.  Contours  in  the  q-plane  for  Im(h1)  =  constant  (from 
Reference  14). 


175 


ESD-TR-81 -1 02 


Appendix  A 


Figure  A-3.  Contours  in  the  q-plane  for  R  (h,)  =  constant 
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f(q)  =  e 


2/3  jq3/2 


( A-4 ) 


2  .  3/2 
-3  3q 


( A-5 ) 


where 


11  1 

„  „3  6  2 

A  =  2  e  it 


=  0.853667218838951 


(9-4)  (81-4)  .  .  .  [9(2  m-1  )  -  4] 

24m  3m  m! 
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Then : 


h  (q)  ~  f  ( q)  u(q),  0  <  arg(q)  <  tt 


( A-6) 


and 


h2(q)  ~  g(q)  v(q)  + 


0  /  0  <  arg(q)  <  ir/2 

e  4tt  j  /  3  f(q)  u(  q)  ,  i  <  arg(q)  <  it 


(A- 7a) 


(A-7b) 


The  reason  that  different  asymptotic  expansions  are  needed  in  different 
regions  of  the  complex  plane  (as  in  Equation  A-7)  is  the  fact  that,  although 
h.j(q)  and  h2(q)  are  analytic  for  all  finite  values  of  q,  these  functions  have 
branch  cuts  at  infinity.  Since  the  asymptotic  expansion  is  essentially  an 
expansion  about  a  point  at  "infinity",  the  branch  behavior  becomes  apparent  in 
the  asymptotic  expansion.  As  discussed  in  Reference  14,  Equation  A- 7a  is 
valid  in  the  region  shown  in  Figure  A-5a,  while  Equation  A- 7b  is  valid  in  the 
region  shown  in  Figure  A-5b.  lb  avoid  approaching  the  branch  cut,  it  is 
convenient  to  use  the  region  division  used  in  Equation  A-7. 

From  Equations  A-1,  A-6  and  A-7: 


F(q)  ~  g(q)  v(q) 


e47Ti/3  f  (q)  u(q)  ,  0  <  arg(q)  <  — 

0  '  \  <  ar<?<!?)  <  n 


( A-8 a) 

( A-8b ) 
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Im(q) 


im(q) 


(a) 


(b) 


Region  of  validity  of 
Equation  A- 7a 


Region  of  validity  of 
Equation  A- 7b 


Figure  A-5.  Regions  of  validity  of  asymptotic  expansions  for  h^. 


Tile  exponential  behavior  of  h ^  and  F  will  enter  through  the  variables 
f  and  g  as  defined  by  Equations  A-4  and  A-5.  Letting: 


j  0 

q  =  p  e 


(A-9) 
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where  p  and  0  are  real,  would  produce: 


3/2  3/2  ,  3  .  .  .  3 

I  =  p  (cos  —8+3  sin  —  0) 


(A-1 0) 


so  that 


_  j  8  j6 

f  =  a  eJ  ,  g  =  y  e 


( A-1 1 ) 


where 


a  =  e 


23/2  .  3  c 

3  p  Sln  7  9 


Y  =  e 


2  3/2  .  3  . 

3  P  Sln  2  9 


.  23/2  3  „ 

0  =  ?  P  cos  2  9 

,  23/2  3  . 

£  =  -  —  p  cos  —  0 


( A-1 2 ) 


and  a,  0,  y  ana  5  are  real  numbers.  Since  the  modulus  of  and  e^  is 
unity,  the  magnitude  of  f  and  g  may  be  obtained  from  a  and  y,  respectively. 
For  a  particular  value  of  p  =  jq|,  a  and  y  a^e  determined  by 
0  =  arg  (q).  Whether  a  and  y  are  large  or  small  would  depend  on  the  sign  of 
the  exponents  in  equation  A-1 2,  which  in  turn  depends  on  the  sign  of 
sin  -j  0.  It  is  clear,  though,  that  when  a  is  exponentially  large,  y  is 
exponentially  small,  and  vice  versa.  Now: 
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sign 


(sin  |  0) 


0  <  8  <  2it/3 
2tt/3  <  8  <  ir 


( A-1  3  ) 


where  only  the  values  of  G  <  0  <  it  are  considered  since  only  the  upper  half  of 
the  q-plane  is  of  interest.  Therefore,  for  p  large: 


a  <<  1 , 

Y  >>  1 ,  0  <  8  <  -^ 

A 

A 

,  2ir 

Y  <<  1  r  —  <  6  <  TT 

or 

l f i  «  i, 

| g |  >>  i ,  o  <  e  <  -^ 

1 f 1  >>  1r 

|g  1  <<  1  r  ~  <  0  <  It 

Therefore,  from 

Equations  A-6,  A-7  and  A-8, 

(A-1 4) 


(A-1 5) 


1  h  1  ( q)  |  <<  1, 

_  ,  ,  2it 

0  <  arg  (q)  <  — 

|h,(q)|  >>  i, 

2  it  , 

—  <  arg  (q)  <  it 

(A-1 6) 


I  h  2<  q) I  >>  Ir 


0  <  arg  (q)  <  it 


(A-1 7 ) 
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| F(q) |  <<  1, 

where,  in  Equations  A-7b  and  A-8a,  the  fact  was  used  that  the  sum  of  an 
exponentially  large  term  and  an  exponentially  small  term  is  expontentially 
large.  The  inequalities  in  the  above  equations  should  be  interpreted  in  terms 
of  the  "potential"  for  being  much  less  than  or  greater  than  1,  since  it  is 
clear  that  as  the  boundary  0  =  arg(q)  +  2tt / 3 ,  the  inequality  would  not 
necessarily  hold. 

Prom  Equations  A-16,  A-1 7  and  A-1 8,  it  is  seen  that  | h 1 (q ) |  and  |h  (q) | 
may  both  be  exponentially  large  simultaneously  in  a  portion  of  the  upper  half 
plane,  whereas  h1 (q)  and  F(q)  cannot  be  large  simultaneously  in  this  region. 
This  fact  prompts  the  use  of  the  functions  (q)  =  h^ (q)  and  K^(q)  =  F(q)  in 
the  solution  to  the  Stokes  equation. 


0  <  arg  (q)  < 


2  ir 

—  <  arg  (q)  <  tt 


Appendix  A 

(A-1  8  ) 
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APPENDIX  B 

EVALUATION  OF  THE  MODAL  DETERMINANT 


GENERAL 

In  searching  for  the  roots  of  the  modal  determinant,  this  determinant 
must  be  evaluated  many  times.  An  efficient  method  for  accomplishing  this  is, 
therefore,  required.  Standard  elimination  methods  require  many  summing 
operations  which  should  be  avoided  as  much  as  possible  if  there  are  large 
differences  in  magnitude  among  the  elements  of  the  determinant.  Ihe 
evaluation  method  should  also  take  maximum  advantage  of  the  presence  of  zero 
elements  in  the  determinant.  A  method  is  described  below  to  accomplish 
this.  This  will  be  done  first  for  a  simple  5-by-5  matrix.  It  will  then  be 
generalized  to  the  N-by-N  case. 

5-by-5  MATRIX 


When  there  are  three  atmospheric  layers  (L  =  3),  the  modal  determinant  is 
given  by  Equation  136  which  is  written  here  as: 


a11 

ai  2 

0 

0 

0 

a21 

a22 

a23 

a24 

0 

a31 

a32 

a33 

a34 

0 

0 

0 

a43 

a44 

a45 

0 

0 

a53 

3  54 

a55 

It  is  desired  to  evaluate  this  determinant. 
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The  determinant  in  Equation  B-1  can  be  written  as  a  linear  sum  of 
cofactors : 


|A|  =  an  M(an)  -  a12  M(a12) 


(B-2) 


where  the  fact  that  a  =  a  =  a  =  0  was  used,  and  the  determinants  of  the 
minors  Mfa.^)  and  M(ai2)  are: 


a22 

a23 

a24 

0 

a32 

a33 

a34 

0 

"<*n>  ’ 

0 

a43 

a44 

a45 

0 

a53 

a54 

a55 

0 

3  22 

3  23 

a24 

0 

9  32 

3  33 

3  34 

in 

0 

a43 

a44 

a55 

0 

a53 

3  54 

( B-3 ) 
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a2l 

a23 

a24 

0 

a21 

0 

3  23 

3  24 

a3l 

a33 

a34 

0 

a31 

0 

3  33 

3  34 

= 

0 

a43 

a44 

a45 

= 

0 

a45 

a43 

a44 

0 

a53 

a54 

a55 

0 

3  55 

a53 

3  55 

(B-4b) 


The  final  determinants  in  Equation  B-4  were  obtained  by  permuting  the 
columns  of  the  respective  matrices. 


Consider  a  system  of  four  equations  in  four  unknowns : 


Bu  =  C 


(B-5) 


where  B  is  the  known,  non-singular  matrix  formed  from  A  as  the  minor  of  the 
a15  term: 


u  is  the  unknown  vector 
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and  C  is  the  known  free  vector  formed  by  the  negative  of  the  last  column  of  A 
excluding  the  first  element: 
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Mta, 1 ) |  =  |M(ai5) |  Ui 


(  B — 1  1  ) 


and 


M(a, 2> |  =  -  |M(a15) |  u2 


( B-1 2) 


so  that  Equation  B-2  may  be  written: 


|M(a15) 


U11U1 


ai2U2) 


(  B-1 3) 


From  Equation  B-6,  it  is  seen  that  the  determinant  of  M(a15)  is  simply 
the  product  of  the  determinant  of  the  2-by-2  submatrix  in  the  upper  left 
corner  and  the  determinant  of  the  2-by-2  submatrix  in  the  lower  right 
corner.  That  is: 


a  ^ 

a 

a  - 

a  „ , 

21 

22 

43 

44 

a„ 

a^ 

• 

a. 

a_ 

31 

32 

53 

54 

Substituting  Equations  B-6,  B-7,  and  B-8  in  Equation  B-5  yields  the 


system  of  equations: 
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21 


+  a„„  u„  +  a„„  u„  +  a_  u , 
22  2  23  3  24  4 


0 


31 


+  a32  U2  +  a33  U3  +  a34  U4 


0 


a  u  -fa,  u 
43  3  44  4 


45 


a  u 
53  3 


a54  U4 


55 


It  is  clear  from  Equations  B-15  that  the  values  of  u3  and  u4  may  be 
by  solving  Equations  B-15c  and  B-15d  alone: 


a4  5 

a44 

-a_  _ 

a 

55 

54 

a43 

a44 

d53 

a54 

a43 

“a45 

a53 

-a55 

a43 

a44 

a53 

a54 

where  it  is  assumed  that  the  denominations  have  non-zero  values. 


( B-1 5a) 


(B-1 5b 1 


(  B-1 5c) 


(  B-1 5d ) 


obtained 


(B-16) 


(B-1 7) 
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Equations  B-16  and  B-17  may  now  be  substituted  into  Equations  B-l 5a  and  B-15b 
to  obtain: 


a  _  u 
21  1 


+  u„ 

22  2 


w 

1 


a43 

a44 

a 

a  „ 

53 

54 

C  B- 1  8 ) 


31  1 


+  a 


32 


43 

a44 

53 

3  54 

(B-19) 


where 


a45 

a44 

+ 

a  ^  , 

3  43 

a45 

a55 

a54 

24 

a53 

a55 

( B-20) 


and 


= 

a45 

3  44 

+  a  ^  , 

a43 

3  45 

33 

a55 

354 

34 

a53 

a55 

Equations  B-18  and  B-19  may  be  easily  solved  to  yield: 


u 


1 


22 
1  32 


|M(ai5) 


( B-22 ) 
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where  Equation  B-14  was  used. 


Substituting  Equations  B-22  and  B-23  in  Equation  B-1 3  results  in: 


w 

a  _ 

a 

w 

1 

22 

21 

1 

+ 

a .  „ 

a  _ 

1 2  a 

w 

2 

32 

31 

2 

where  w1  and  w^  are  given  by  Equations  B-20  and  B-21,  respectively. 

N-by-N  MATRIX 

The  method  utilized  above  for  a  5-by-5  matrix  will  now  be  generalized. 
Consider  the  N-by-N  matrix  given  by: 
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From  the  definition  of  the  determinant  of  a  matrix  in  terms  of  a  cofactor 
expansion : 


Z  (-1  )  a  |M  (a  ) I 
i=1  1  i  1i 


( B-26 ) 


Define  the  following  permutation  operation  on  the  columns  of  a  square 
L 

L-by-L  matrix  E:  Let  P^(E)  be  the  matrix  formed  by  the  operation  of  moving 
columns  a,  a+1 ,  a+2,  .  .  .  ,L-1,  L  so  that  column  L  becomes  column  a,  column 
a  becomes  column  a+1,  column  a+1  becomes  column  a+2,  etc.  Now: 


PL  (E) 
a 


(-1  ) 


L-a 


( B-27  ) 


In  particular,  let 


E  =  M( a  .  ) 

1i 


( B-28 ) 


so  that  L  =  N  -  1,  and  let  a  =  i.  Thus: 


=  (-1  ) 


N-1-i 


P. 

i 


N-1 


(E)  | 


( B-29  ) 


Then 


N-1  i+1  N-1 -i  N-1  ,  N+1 

A  =  Z  (-1)  a  (-1)  P  (M( a  )>  +  (-1)  a  |M(a  )| 

i  =  1  1i  i  1  i  IN  IN 

( B-30 ) 

N  N-1  ,  N-1  ,  N+1 

=  (-1  )  Z  a  |  P  (M(a  ))  |  +  (-1  )  a  |  M( a  )  | 

i= 1  1i  i  1 i  IN  IN 
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Now  consider  the  system  ot  N-t  equations  and  N-1  unknowns  : 


Bu  =  O  ( B- 3 1  ) 


where  B  is  the  known  matrix  formed  f  r  cm  A  as  the  minor  ot  a,.,: 

1  N 


B 


M<  a  ) 
IN 


'’n,  N-1 

(B-JjT 


u  is  the  unknown  vector 


( B-333 


and  r  is  the  known  trot'  vector  formed  by  the  negative  ot  the  last  column  of  A 
excluding  the  first  element: 


l».j 


( B-34) 
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It  is  known  that,  if  the  matrix  B  is  non-singular,  the  solution  for  u^  is 


obtained  as: 


IB^I 


,  1  <  i  <  N-1 


(B-35) 


so  that 


|B(i)|  =  u.  | B |  , 


1  <  i  <  N-1 


( B-36) 


where  B^)  is  the  matrix  obtained  by  replacing  the  ith  column  of  B  with  the 
free  vector  C.  But  some  thought  would  show  that: 


|B(i)|  =  -  IP.**"1  (M  (au))|. 


1  <  i  <  N-1 


( B-37 ) 


Using  Equations  B-32  and  B-37  in  Equation  B-36  leads  to: 


|P.  ( M ( a  ) ) I  =  -  u.  | M{ a  ) | ,  1  <  i  <  N-1  (B-38) 

1  11  1  IN  ~  ~ 


Substituting  Equation  B-38  into  Equation  B-30  yields: 

N  N-1  N+1 

I A  |  =  (-1)  {-  l  a  u  |M(a  )|}  +  (-1)  a  |M(a  )| 

i=1  1i  i  IN  IN  IN 


N+1  N-1 

=  (-1 )  | M(a  )  {a  +  l  a  u  } 

IN  IN  i=1  1i  i 
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or 


N+1 

(-1)  S  M(  a  )j 
IN 


N 

la  u  , 
i=1  1  i  i 


u 

N 


( B-39 ) 


APPLICATION 

At  first  glance,  Equation  B-39  looks  more  involved  than  its  equivalents 
in  Equations  B-26  and  B-30.  Indeed,  the  form  of  Equation  B-39  contains  the 
unknowns  in  the  vector  u,  which  first  must  be  solved  through  Equation  B-31 
before  a  solution  for  } A }  is  obtained.  Nevertheless,  there  exist  matrices  for 
which  Equation  B-39  represents  a  more  efficient  method  for  evaluating  )a|  than 
does  Equation  B-26.  Such  a  matrix  would  be  one  for  which  B  =  M(  a  )  would 
be  relatively  simple  to  evaluate  even  though  ]a)  is  not. 

Consider,  for  example,  the  N-by-N  matrix  A  for  which  the  submatrix  B  = 

M(a,N)  has  the  form  illustrated  in  Figure  B-1.  About  the  diagonal 

of  M(a,„)  are  2-by-2  and  3-by-3  submatrices  which  are  each  (except  for  those 
IN 

on  the  array  boundary)  flanked  below  and  to  the  left  by  zero  elements.  There 
are,  say,  K  such  submatrices.  The  kth  such  matrix  will  have  N^-by-N^ 
elements,  so  that: 

K 

L  N  =  N-1  ( B-40 ) 

k=1  k 


Define  the  determinant  of  as: 


\  ■ 


(B-41) 
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The  rectangular  matrix  consisting  of  the  elements  of  B  to  the  right  of  will 
be  called  the  matrix  F^,k  =  1,2  K-1.  The  elements  of  each  F^  may  be 

zero  or  non-zero.  F^  will  have  N^  rows  and  N^_  columns,  where: 


-  K 

N  =  E  N  ( B-42 ) 

k  i=k  +  l  i 


/  k  i 

Now  define  the  subvectors  u'  '  of  the  unknown  vector  u  as  that  portion  of 

( k  1 

u  corresponding  to  the  matrix  F^  (see  Figure  B-1 ) .  uv  '  will  contain  N^ 

terms.  Its  first  term  will  be  the  n^1*)  element  of  u,  and  its  last  term  will 

/ 1.  \ 

be  element  number  np'  '  of  u,  where 


k-1 

(k)  E  N  ,  k  >  1 

n  =14-  i=i  i  ( B-43 ) 

I 

0  ,  k  =  1 


and 


(k)  k 

n  =  E  N 

F  i  =  1  i 


( B-44 ) 


I  k  1 

Similarly,  define  the  subvectors  C'  '  of  the  free  vector  C. 

_  ()t ) 

Also  define  complementary  subvectors  u  of  the  unknown  vector  u. 

—  ( k) 

u  is  the  vector  consisting  of  all  terms  in  u  beyond  those  which  make  up 
u<k>.  Thus  the  first  term  in  u  will  be  the  term  number  n  of  u  and 

the  last  term  will  be  term  number  N-1  of  u,  where: 
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no 


+  1 


( B-45 ) 


It  is  now  clear  that: 


|M  (a  ) 
IN 


K 

n  d 

k 


( B-46 ) 


( K ) 

Now  the  subvector  ul  may  be  solved  using  the  equation: 


«  u™  » 
K 


(B-47) 


so  that 


u. 

i 


<K) 


( B-48 ) 


where,  again,  M  is  the  matrix  obtained  by  replacing  the  ith  column  of  Mv 

K.  In 

with  the  vector 

( K— 1  ) 

The  subvector  u  may  be  solved  from  the  equation: 


.(K-1 ) 


K-1 


C(K-1) 

C  (K'1) 


F  u(K_,) 
K-1 


( B-49 ) 


so  tha  t 
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u 

l 


(K-1  ) 


(i) 


K-1 


K-1 


( B-50) 


where  M  (i)  is  the  matrix  obtained  by  replacing  the  ith  column  of  M  by 

K— '  _ (k  i )  K —  1 

the  vector  C  .  The  kth  subvector  u^  '  may  now  be  solved  for  through  the 

equation : 


\ 


:(K)  =  C(k)  -  F  u  (k),  c  (k) 

k 


so  that 


(B-51  ) 


u. 

i 


(k) 


(B-52) 


for  the  subvector 
Equation  B-46  in 

SINGULAR  SUBMATRIX 


This  procedure  may  be  continued  until  a  solution  is  obtained 
u*1^.  The  values  of  u^  so  calculated  may  be  used  along  with 
Equation  B-39  to  obtain  the  determinant  of  A. 


The  effect  of  one  of  the  submatrices  (say  )  being  singular  is  now 
considered.  In  solving  this,  an  investigation  is  made  into  the  manner  in 
which  Dj  enters  the  final  expression  for  A  when  is  small. 

It  is  clear  that  would  not  enter  the  successive  calculations  described 
in  Equations  B-51  and  B-52  until  the  value  of  u^ ^  is  to  be  evaluated.  Then: 


u. 


(j) 


D. 

1 


D, 

3 


( B-53 ) 
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where  the  vector  '  is  defined: 

/ |Mi  '\ 

v(^  = 

/  '  \ 

|m. (2) |  \ 

(  B-5 

|mJ(3)  | 

l 

t 

Vov 

<* 

The  equation  for  ^  is  then: 
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B 


Substituting  this  into  Equation  B-55  yields: 


M.  u 
1-1 


j-1 


1  (D  C(j_1)  -  F  v 


°1 


j-1 


(B-57) 


As  Dj  becomes  small,  the  first  term  in  the  brackets  of  Equation  B-57  will 
vanish,  as  will  all  elements  in  the  product  F.  v  ^  1 ^  containing  the  factor 
Dj  .  But  this  is  the  result  that  would  be  obtained  if : 


=  (j+1) 


,(j) 


U 

D. 

1 


(3+2) 


(K) 


(B-58) 


E^ch  term  in  the  free  vector  in  Equation  B-57  will  therefore  be  inversely 

proportional  to  D-,  In  solving  Equation  B-57,  each  element  of  the  vector 

u(j  1 ^  will  likewise  be  inversely  proportional  to  .  The  same  reasoning  may 

be  used  for  u^  ^ u^1^  to  show  that  each  of  these  subvectors 

will  be  inversely  proportional  to  D- ,  and  that  the  result  obtained  requires 

(k )  ( ^ ) 

ignoring  the  subvectors  C  and  u  ,  j  +  1  j<  k  _<  K,  in  comparison  to  the  other 
terms  of  the  free  vector.  Using  Equations  B-58  and  B-46  in  Equation  B-39 
yields : 


N+1 

(-1  ) 


K  N 

HD  £  a  u 

k=1  k.  i=1  1i  i 


N+1 

=  (-1  ) 


K 

n 

k=1 


k^j 


D 

k 


N 

E  a  u 
i=1  1i  i 


(B-59) 
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where 


( B-60 ) 


in  which  the  factor  in  the  numerator  and  the  denominator  have  been 
cancelled,  and  does  not  appear  in  Equation  B-59.  The  unknown 
u(i),  1  £  i  _<  j-1  are  assumed  to  be  obtained  using  the  successive  method 
described  by  Equations  27  and  28,  using: 


u 


0 


0, 


j  +  1  <  k  <  K 


(  B-61  ) 


and  with  u^  ^  replaced  by  v^  ^  .  Equation  B-59  is  the  final  value  for  the 
determinant  when  the  j th  submatrix  ft  is  singular. 

In  the  event  that  both  M,  and  M.  are  singular,  where  a  is  a  positive 

3  3  -a 

integer,  then  the  resulting  value  of  | A |  would  be: 


N+1  K 

(-1  )  n 

k  =  1 


N 

D  Z  a  u 
k  i=1  1i  i 


k* j ,  j -a 


( B-62 ) 


where 
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where 


(  k) 

and  the  u  ,  1  <  k  <  j-  a  -  1  are  obtained  as  described  above.  The 
extension  to  any  number  of  singular  submatrices  is  straightforward. 


PARTICULAR  CASE 


Of  particular  interest  will  be  the  matrix  A  with: 


a  =  0,  N  +  1  £  i  _<  N 


( B-64 ) 


a .  =0,  N-N  +  1  <  i  <  N  (B-65) 

IN  K  —  — 


and  with 


(Fk)  =  °' 

a6 


8  > 


k+ 1 


( B-66 ) 
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N 

N  +  1  1 

A|  =  (-1)  | M( a  )|  l  a  u 

IN  i=1  11  i 


(B-67) 


Prom  Equation  B-65: 


C 


(  1  ) 


,  ( 2 ) 


C(K-1) 


0 


(  B-68 ) 


and  from  Equation  B-66,  the  submatrix  M(  a  )  is  as  shown  in  Figure  B-2.  The 

1  N 

non-zero  portion  of  the  F^  will  be  called  and  will  have  rows  and 

columns. 

The  first  step  in  the  solution  procedure  is  determination  of  the  elements 
of  u'  '  from  Equations  B-47  and  B-48  to  yield: 


u 


(K) 


{ B-69 ) 


(  K  ) 

where  again  v'  is  the  vector,  the  ith  term  of  which  is  the  determinant  of 

( K ) 

the  matrix  formed  by  replacing  the  ith  column  of  by  the  free  vector  C'  . 
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From  Equations  B-49  and  B-68,  the  elements  of  u^-1  ^  are  given  by: 


M  u 
K-1 


(K_1)  -  -F  u^-15 

K-1 


fk-i  u 


(K) 


(K) 


=  F 


K-1  D 


c  {K-'] 


(B-70) 


This  has  the  solution: 


u 


(K-1  ) 


(K-1  ) 


D  D 
K  K-1 


(  B-71  ) 


( K  —  1  ) 

where  now  v'  is  the  vector,  the  ith  term  of  which  is  the  determinant  of 

_  (K-1 ) 

the  matrix  formed  by  replacing  the  ith  column  of  MK_1  by  the  vector  C 

f  k  1 

The  elements  of  u'  '  are  given  by: 


(k)  •  (k) 

Mu  =  -F  u 
k  k 


(k-1  ) 

’  v 

F  - 

k  00  0  «  •  •  D 

K  K-1  K-2  k-1 


=  -F 


(k-1) 


K 

n  d 

i  =  k-1  i 


-  (k) 
C 

K 

n 

=  k-1 


( B-72 ) 


with  the  solution: 
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(k) 

u 


(k) 

v 


K 

II  D 
i  =  k  i 


lhe  solution  of  1  ^  will  similarly  be  given  by: 


( B-73) 


(1  ) 

(1  )  v 

u  =  - 

K 

n  d 
i  =  1  i 


Substituting  Bguation  B-46  and  B-74  into  Equation  B-67  yields: 


( B-74 ) 


N 

.  ,  N+1  1  (1) 

A  =  (-1)  Z  a  v 

i=1  1 i  i 


( B-75 ) 


Thus  it  is  seen  that,  for  this  particular  form  of  the  matrix  A,  the 
determinants  of  the  submatrices  do  not  enter  into  the  final  expression  for 
| A | .  It  follows  that  Equation  B-75  also  holds  when  any  of  the  submatrices  is 
singular . 

EXAMPLE  -  FIELDS  IN  HOMOGENEOUS  STRATIFIED  LAYERS 


Consider  L  layers  of  atmosphere  over  a  ground,  with  each  layer 
homogeneous  and  characterized  by  a  parameter  a.  In  acoustics,  this  parameter 
might  be  the  medium  density.  In  electromagnetics,  this  might  be  the 
permittivity  e.  A  field  quantity  E^  in  the  ith  layer  may  be  written  as: 


E^  (r,  z)  =  /  f  (cr  k,  z)  f  (k,r)  dk 


( B-76 ) 
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where  f  is  a  known  kernel  and  the  integration  is  over  a  closed  contour  in  the 
complex  k-plane.  Assuming  the  quantity  ¥  satisfies  a  linear  second-order 
differential  equation,  it  may  be  written  as: 


(a. ,k,z)  =  A.  h  ( a . , k , z )  +  B.  h  (a.,k,z)  +  <5.  g(a. ,k,z,z  )  (B-77) 

ii  ill  i  2  i  ipi  T 


where 


h^  and  h^  are  linearly  independent  functions, 

g  is  a  Green  function,  which  is  a  particular  solution  of  the  governing 
differential  equation  that  accounts  for  a  source  at  the  location 
z  =  z^  located  in  the  pth  layer, 

ip  is  the  Kronecker  delta  function. 


The  unknowns  A^  and  B^  may  be  found  by  solving  the  system  of  linear  equations 

representing  the  boundary  conditions  on  the  ?.  and  f , '  =  dY./dz  at  the  layer 

i ,  i  i 

interfaces  (z=z,1<i<L-1)  and  at  the  ground,  (  z  =  z  )  and  the 
i  —  —  o 

radiation  conditions  at  z  ♦  +  •.  For  L  =  4  This  system  will  have  the  form: 


A1  Vq10>  +  B1  h2<<ll0) 

A1  h1^11^  +  B1  h2(qi  1 }  ~  A2  “  B2  h2^21^ 

Ai  hi  +  Bi  h2^11^  “  a2  l,l(£l21^  "  B2  ^2^21^ 

A2  hi(q22)  +  ^2^22^  “  A3  ^1^32^  “  B3  ^2^32^ 

A2  h{(<l22^  +  b2  h2^22^  “  a3  hl^32^  “  B3  h2^q32^ 

A^  h^j  (qjj)  +  B3  ^2 ( q3 3 )  “  B4  ^2^^43^ 

a3  hl(<l23)  +  b3  h2^33^  “  b4  h2^43* 


(B-78) 


where  h  (q . . )  =  h  (a.,k,z),  j*i,  i+1 ,  i  >  0,  h  (q  )  =  h'(q  )  -  Gh  (q  ) , 
m  ji  m  1  j  m  TO  m  (0  m  10 

where  G  is  a  function  of  k  and  the  ground  parameters,  and  it  has  been  assumed 


that  p=2.  The  f^  are  functions  of  k.  Equation  B-76  may  be  evaluated  using 
residue  theory,  where  the  poles  of  are  the  values  of  k  at  which  the 
determinant  of  the  above  system  of  equations  vanishes.  When  the  roots  of  this 


209 


ESD-TR-81-1 02  Appendix  B 

determinant  are  evaluated  numerically,  the  determinant  must  be  evaluated  many 
times,  thus  making  an  efficient  means  for  such  evaluation  very  useful. 

The  matrix  has  the  form: 


This  matrix  has  the  form  of  the  matrix  in  Figure  B-2,  and  the  evaluation  of 
its  determinant  may  be  accomplished  using  the  algorithm  described  above.  The 
total  number  of  multiplication  operations  required  is  on  the  order  of  8K  where 
K  is  the  number  of  2nd-order  submatrices.  In  the  example  above,  K=  3,  so  that 
the  number  of  multiplication  operations  is  24,  or  about  the  number  of  non-zero 
terms  in  the  matrix.  This  compares  with  N!  for  the  case  of  expansion  of  the 
determinant  in  cofactors,  which,  for  the  above  example,  would  be  over  5000. 

In  the  event  the  submatrices  were  of  order  3,  then  the  number  of 
multiplication  operations  would  be  on  the  order  of  45K.  If  K=  2,  in  a  7-by-7 
matrix  such  as  Equation  B-79,  this  number  would  be  90,  which  still  is  small 
compared  with  the  5000  operations  required  to  evaluate  the  determinant  of  a 
7-by-7  matrix  using  a  cofactor  expansion. 
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